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The  Optimal  Design  of  Structures  Incorporating  ‘Smart  Materials’ 


Current  Objectives 


1 .0  Current  Objectives 

The  purpose  of  this  project  was  to  develop  a  rigorous  theoretical  foundation  of  the  design  of 
structures  incorporating  general  ‘smart'  materials,  and  to  explore,  extend  and  develope  the 
use  of  sensors  and  actuators  that  employ  ‘finite  length’  paths.  It  is  anticipated  that  this 
project  will  fundamentally  impact  the  practice  of  designing  passive  and  active  materials,  , 
and  structural  topologies.  \ 

There  were  three  main  tasks  to  be  performed  under  this  grant.  These  tasks  are  listed  below 
along  with  some  of  the  details  of  our  current  approach  and  status.  Even  though  this  grant 
has  ended,  some  of  this  work  (as  described  below)  is  continuing  with  support  from  various^ 
other  sources. 

These  objectives  are  identical  to  those  that  were  listed  in  the  proposal  and  the  last  progress 
report  with  one  exception:  an  additional  task,  ‘practical  implementation  of  the  optimization 
procedure  for  large  scale  problems’  has  been  appended.  As  will  be  briefly  discussed  below 
and  borne  out  through  most  of  the  references,  in  some  cases  the  objectives  have  been 
achieved  and  in  others  we  have  at  least  made  substantial  progress.  That  is,  at  a  minimum  we 
have  been  able  to  create  simple  problem  formulations  that  take  into  account  each  of  the 
tasks  (e.g.  General  Design  Objectives,  General  Local  Design,  General  Tensor  Design,  ect.). 

1.1  Optimal  Design 

The  purpose  of  this  main  task  is  to  develop  a  rigorous  theoretical  foundation  to  design  structures  incorporating 
‘smart’  materials.  The  aspects  to  be  explored  include: 

1.1.1  General  Design  Objectives 

Rather  than  be  restricted  to  designing  for  ‘minimum  compliance’,  other  objectives  such  as  ‘maximum 
strength1,  or  ‘maximum  life’  within  the  context  of  a  combined  loading  environment  will  be  examined. 

1.1.2  General  Local  Design  Constraint 

The  current  formulation  has  very  limited  means  to  restrict  local  type  of  material  being  optimized.  Here,  a 
formulation  that  allows  the  local  material  to  be  restricted  in  some  manner  (e.g.  isotropic  or  to  be  ‘man¬ 
ufacturable’)  will  be  developed. 

1 .1 .3  General  Tensor  Fields 

The  present  analysis  ‘designs’  the  constitutive  tensor  of  a  generally  linear  anisotropic  material.  Here  oth¬ 
er  tensor  fields,  such  as  those  associated  with  ‘smart’  materials,  will  be  incorporated  into  the  design. 

1.1.4  Nonlinear  Materials 

Many  of  the  above  elements  to  be  incorporated  will  lead  to  certain  non-linearities  in  the  design  process. 
Here  the  aim  is  to  investigate  the  design  of  materials  that  have  constitutive  non-linearities. 

1.1.5  Nonlinear  Kinematics 

Many  current  ‘solutions'  to  aerospace  problems  that  employ  smart  materials,  are  mechanisms.  It  is  pos¬ 
sible  to  incorporate  nonlinear  kinematics  and  potentially  ‘contact’  mechanics  into  the  optimal  design  for¬ 
mulation. 
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1.1.6  Practical  Implementation 

How  the  numerical  implementation  of  these  optimization  routines  scale  to  larger  more  practical  problems 
is  not  obvious.  In  fact  it  appears  that  to  solve  problems  of  any  type  of  substantial  size  may  require  alot  of 
computation  effort  (e.g.  parallel  processing  machines)  or  some  imaginative  solution  procedures.  This  is 
being  listed  as  a  separate  task  because  this  difficulty  spans  all  the  above  optimization  tasks. 


1.2  Distributed  Sensing  and  Actuation 

The  goal  here  is  to  further  extend  and  explore  aspects  of  sensors  and  actuators  that  employ  ‘finite-length’ 
paths. 


1.2.1  Theoretical  Investigation 

The  current  isotropic  analysis  is  to  be  used  to  explore  bounds  on  the  feasibility  of  using  these  types  of 
sensors.  The  analysis  is  to  be  extended  to  include  an  explicit  measure  of  structural  integrity  sensing.  It  is 
also  to  be  extended  to  more  general  materials  (e.g.  the  anisotropic  material  of  a  composite  torque  tube) 
and  configurations  (e.g.  a  rotor  blade). 

1 .2.2  Prototype  Construction  and  Testing 

In  cooperation  with  Boeing  Helicopter,  a  prototype  torque  tube  is  to  be  constructed  and  bench  tested.  The 
ability  of  the  system  to  measure  torque,  and  structural  integrity  will  be  examined. 

1.2.3  Practical  Implementation 

The  current  means  to  ‘design’  distributed  sensors  involves  knowing  an  analytical  model  of  the  system.  A 
goal  of  this  aspect  of  the  investigation  is  to  ‘design’  distributed  sensors  and  actuators  (using  the  approach 
of  the  ‘optimal  design’  above)  on  more  practical  structures. 

1 .2.4  Performance  of  an  Electroactive  Hydraulic  Actuator 

Certain  electroactive  materials  (Re.  Shoko  Yoshikawa)  show  a  substantial  volume  change  under  an  ap¬ 
plied  magnetic  field.  Under  certain  limiting  conditions  similar  to  those  proposed  by  other  members  of  the 
Smart  Structures  Rotorcraft  Consortium  these  materials  can  potentially  achieve  a  1%  change  in  volume. 
This  is  an  effect  that  is  huge  (almost  an  order  of  magnitude)  greater  than  the  intrinsic  effects  being  pro¬ 
posed  for  other  actuator  designs.  The  purpose  here  it  to  construct  a  small  prototype  pump  that  will  allow 
feasibility  of  using  the  volumetric  effect  of  certain  electroactive  materials  in  hydraulic  systems. 


1.3  Interaction  with  Other  Program  Participants 

We  will  cooperate  with  other  participants  so  as  to  be  in  a  position  to  disseminate  results,  incorporate  their  ob¬ 
jectives  and  constraints  into  this  study. 
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2.0  Status  of  Effort 


Substantial  progress  has  again  been  made  toward  fulfilling  all  the  objectives.  As  alluded  to 
in  the  last  report  we  have  extremum  principles  that  now  work  for  all  of  the  tasks  listed 
above.  Listed  below  is  a  description  of  the  status  in  each  of  the  areas.  In  some  cases  there  is 
reference  to  publish  documents  or  the  papers  appended  at  the  end  of  this  report. 

2.1  Optimal  Design 

In  general  we  have  addressed  or  solved  at  some  level  every  topic  here.  However,  in  some  cases  these  issues 
will  remain  as  research  topics  for  the  foreseeable  future.  The  topics  being  explored  or  that  have  been  solved 
include: 

2.1.1  General  Design  Objectives 

Designing  for  ‘minimum  compliance’  is  still  used  as  a  baseline  objective.  However  other  objectives,  es¬ 
pecially  at  the  urging  of  some  of  the  industrial  contacts  (e.g.  ALCOA),  such  as  ‘minimum  cost’  or  ‘min¬ 
imum  life-cycle  cost’  are  also  being  explored.  To  see  example  of  these  please  refer  to  the  documents  [1, 
5,  61.  There  remains  much  work  to  be  done  here.  For  instance  the  ‘maximum  strength’  objective  has  been 
explored  -  some  simple  sample  problems  (and  still  yet  to  be  published)  have  been  solved. 

2.1.2  General  Local  Design  Constraint 

The  original  formulation  had  very  limited  means  to  restrict  the  local  type  of  material  being  optimized. 
The  original  formulation  used  purely  4-tensor  invariants.  Thus,  any  arbitrarily  rotation  of  the  material  co¬ 
ordinate  frame  could  never  result  in  a  cost.  While  this  approach  has  appeal  from  an  academic  perspective, 
it  fundamentally  ignores  basic  manufacturing  issues.  In  addition,  certain  manufacturing  costs  that  involve 
orientation  dependance  can  also  be  incorporated.  We  can  now  restrict  the  material  to  be  anything  from 
isotropic  to  generally  anisotropic.  Please  refer  to  the  document  [6,  8,  12].  We  can  now  local  restrict  ma¬ 
terial  properties  in  multiple  ways.  This  task  is  essentially  solved  for  our  current  purposes. 

2.1 .3  General  Tensor  Fields 

Optimization  formulations  that  include  distributed  temperature  and  electrical  fields  have  been  formu¬ 
lated.  Here  we  have  a  formulation  that  incorporates  an  elastic  field  and  another  field  that  is  linearly  cou¬ 
pled.  The  interesting  feature  here  is  that  the  solution  drastically  changes  character.  For  instance  for  a 
simple  bar  rather  than  having  simple  ‘diffusion  equations’  (and  for  example  simple  polynomial  solutions) 
you  get  equations  that  are  analogous  to  ‘wave  equations’  (and  for  example  transcendental  solutions).  This 
represent  the  interplay  or  exchange  of  energy  between  the  two  fields.  We  have  some  simple  analytical  so¬ 
lutions  here,  as  well  has  an  extremum  principle!  However,  we  are  currently  working  on  designing  pa¬ 
rameters  (e.g.  the  mixture  of  passive  and  active  materials)  in  this  two  field  system.  Much  of  this  these 
simple  examples  have  yet  to  be  published  however,  the  method  to  separate  the  ‘active’  and  ‘passive’  phas¬ 
es  appears  in  [14,  16]. 

2.1 .4  Nonlinear  Materials 

It  is  a  relatively  straightforward  matter  to  incorporate  nonlinear  constitutive  properties  into  the  typifi- 
cation  formulation.  In  fact,  it  is  perfectly  allowable  to  design  the  nonlinear  portions  of  the  constitutive 
equation.  This  task  is  essentially  completed  for  lots  of  examples.  Please  refer  to  the  documents  [2,  3,  6,  9, 
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2.1.5  Nonlinear  Kinematics 

In  much  the  same  manner  that  nonlinear  constitutive  properties  can  be  incorporated,  nonlinear  kine¬ 
matics  can  also  be  used.  This  task  is  essentially  complete.  Please  refer  to  the  documents  [9,  11] 


2.2  Distributed  Sensing  and  Actuation 

We’ve  essentially  completed  all  that  can  be  done  here  except  the  last,  self-imposed  task. 

2.2.1  Theoretical  Investigation 

A  structural  integrity  sensor  was  analyzed  and  shown  to  exist  for  a  prismatic  homogeneous  isotropic 
structure  with  arbitrary  cross-section.  Current  efforts  are  still  to  extend  this  analysis  to  a  structure  that 
uses  an  anisotropic  material.  Refer  to  document  [4,  7] 

2.2.2  Prototype  Construction  and  Testing 

We’ve  constructed  several  prototypes  which  substantiate  this  approach.  That  is,  it  works!.  Refer  to  doc¬ 
ument  [10]. 

2.2.3  Practical  Implementation 

Since  the  approach  is  technically  sound  we  are  trying  to  locate  an  application  that  requires  the  increase  in 
range  and  stiffness  of  the  transducer.  There  have  been  no  real  buyers  here.  We  have  written  a  proposal 
NASA-JPL  to  make  use  of  this  technique  to  measure  boom  deformations/loads  in  a  space  interferometer 
experiment. 

2.2.4  Performance  of  an  Electroactive  Hydraulic  Actuator 

(This  task  was  not  actually  in  the  original  proposal) 

We’ve  purchased  the  necessary  equipment  to  test  this  device.  However,  we  have  not  had  the  resources  to 
finish  this  task. 


2.3  Interaction  with  Other  Program  Participants 

During  the  course  of  this  grant  there  were  discussions  held  with  many  organizations. 
These  included: 

-  SSRC:  (Boeing/MIT/Penn  State  Consortium) 

-  Aluminum  Company  of  America  (Alcoa) 

-  Ford  Scientific  Research  (Dearborn) 

-  Michigan  Scientific  Company: 

-  Smiths  Industries 

In  addition  numerous  government  agencies  were  contacted  along  with  many  discussions  at 
PI  meetings  at  conferences.  The  most  notable  result  of  all  this  is  that  we  now  have  an 
active  program  with  Ford  Research  and  they  have  incorporated  some  of  our  work  into  their 
research  and  production  design  codes  (e.g.  see  [16]). 
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Also,  in  an  attempt  to  interact  with  the  helicopter  community  we’ve  specifically  addressed 
design  problems  that  involve  inertial  loads  (e.g.  see  [13,15]).  We  anticipate  that  these  will 
be  of  interest  to  Boeing  Helicopter. 

3.0  Accomplishments/New  Findings 


During  this  period  there  were  several  significant  accomplishments.  They  are: 

1 )  We’ve  shown  how  to  include  a  body  force  (for  the  specific  helicopter  rotor  problem)  into 
the  formulation,  and  that  significant  applications  can  be  solved  [13,  15] 

2)  We  can  now  incorporate  an  arbitrary  number  of  materials  (either  passive  or  active)  into 
the  formulation.  Refer  to  documents  [14,  16]. 

3)  A  three-dimensional  version  of  the  optimal  design  code  has  been  implemented  and  ap¬ 
pears  to  be  sufficient  to  handle  problems  of  significant  scale.  Refer  to  document  [15]. 

4.0  Personnel  most  recently  Supported 

Professor  Peter  D.  Washabaugh  ( 1  Summer  Month) 

Professor  John  E.  Taylor  ( 1  Summer  Month) 

Dr.  Sergio  Turteltaub  (9  Months  -  Post  Doc) 

5.0  Past  Report  Publications 


1)  “Analysis  and  Design  of  Trussed  Structures  Made  of  Elastic/Stiffening  Materi¬ 

als”,  J.E.  Taylor  and  P.D.  Washabaugh,  Structural  Optimization,  8  1-8,  1994 

2)  “Applications  of  a  Generalized  Complementary  Energy  Principle  for  the  Equi¬ 

librium  Analysis  of  Softening  Material”,  Plaxton,  S.,  Taylor,  J.E.,  Comp. 
Meth.  Appl.  Mech.  Engng,  1 17,  pp  91-103,  1994 

3)  “A  Global  Extremum  Principle  in  Mixed  Form  for  Equilibrium  Analysis  with 

Elastic/Stiffening  Materials:  A  Generalized  Minimum  Potential  Energy 
Principle”,  Taylor,  J.E.,  Journal  of  Applied  Mechanics,  Vol  61 -No.  4,  pp 
914-918. 

4)  “Elastic  Transducers  Incorporating  Finite  Length  Optical  Paths”,  K.  J.  Peters  and 

P.D.  Washabaugh,  Applied  Optics,  3,  No.  22,  4993-5002,  1995 

5)  “Optimal  Design  of  Material  Properties  and  Material  Distribution  for  Multiple 

Loading  Conditions”,  Bendsoe,  M.P.,  Diaz,  A.,  Lipton,  R.,  Taylor,  J.E.,  hit. 
J.  Num.  Meth.  Engng.  1996 
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6)  “GeneralizedT’otentiar  Formulations  for  Kinematically  'Linear  Elastostaiics  of 

Constitutively  Nonlinear  Structures,  Taylor,  J.E.,  hit.  J.  Solids  and  Struc¬ 
tures,  1996. 

7)  “A  Technique  for  Monitoring  In-Situ  the  Structural  Integrity  of  Prismatic  Struc¬ 

tures  Under  General  Loading”,  K  J.  Peters  and  P.D.  Washabaugh,  AIAA 
Journal,  Vol  35,  No.  4,  1997 

8)  “On  Structural  optimization  Formulations  with  Generalized  Cost  Constraints”, 

Taylor,  J.E.  and  Washabaugh,  P.D.,  Proceedings  PACAM  V Meeting  Puerto 
Rico,  1996. 

9)  “Generalized  Potential  Formulations  for  Elastostatics  of  Constitutively  Nonlin¬ 

ear  Structures  Under  General  Loading”,  Taylor,  J.E.,  International  Journal 
of  Mechanical  Sciences ,  Vol  61 -No.  4,  pp  914-918. 

10)  ‘Application  of  Finite-Length  Displacement  Sensors  to  Precision  Torque  Mea¬ 

surements  of  a  Tail  Rotor  Transmission  Tube”  K.J.  Peters  and  P.D.  Washa¬ 
baugh,  Proceedings  of  SP1E  Smart  Structures  and  Materials  meeting,  San 
Diego  Ca,  1996 

1 1)  “Finite  Strain  Elastostatic  with  Stiffening  Materials  a  Constrained  Minimiza¬ 

tion  Model”,  Hollister,  S.J.,  Taylor,  J.E.,  Washabaugh,  P.  D.,  Journal  of  Ap¬ 
plied  Mechanics,  Vol  63,  No.  2,  1997 

1 2)  “A  Formulation  for  Generalized  Hyperelasticity  using  a  Composition  of  Poten¬ 

tials”,  Guedes,  J.M.,  Taylor,  J.E.,  to  appear  Journal  of  Applied  Mechanics. 
1997 

6.0  Most  Recent  Publications 


13)  “Optimal  distribution  of  material  properties  for  an  elastic  continuum  with  struc¬ 

ture-dependent  body  force”,  Sergio  Turteltaub  and  Peter  Washabaugh.,  To 
appear  in  the  International  Journal  of  Solids  and  Structures,  1998. 

14)  “Prediction  of  Optimal  Material  Layout  (Topology)  and  Properties  for  an  Elas¬ 

tic  Continuum  Structure  using  Weighted  Resource  Constraint”,  J.  Pawlicki, 
J.E.  Taylor,  S.  Turteltaub  and  P.D.  Washabaugh.,  Submitted  to  Structural 
Optimization ,  1998. 

15)  “Three-dimensional  optimal  design  of  structural  topology,  shape  and  material 

properties  with  design-dependent  field  force”,  Sergio  Turteltaub,  Jakub 
Pawlicki,  John  E.  Taylor,  and  Peter  Washabaugh.,  Submitted  to  the  Interna¬ 
tional  Journal  of  Solids  and  Structures,  1998. 

16)  “A  Design  Model  to  Predict  Optimal  Two-Material  Composite  Structures”,  H. 

Rodrigues,  Ciro  A.  Soto.  and.  John  E.  Taylor,  in  preparation,  1998. 


The  current  versions  of  these  papers  are  appended  at  the  end  of  this  report. 
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7.0  Interactions/Transitions 


Please  refer  to  the  section  2.3  for  interactions  with  other  program  participants  and  the  above 
publication  list.  The  most  notable  transition  of  this  information  to  industry  has  come  from 
an  unexpected  source.  In  the  past  year  and  a  half  we  have  been  collaborating  with  Ford  Re¬ 
search  Laboratories.  As  evident  by  the  last  publication  above  [16],  the  work  on  Optimal  To¬ 
pology  Design  has  made  it  into  some  of  the  design  codes  at  Ford.  They  are  continuing  to 
support  this  work. 


8.0  New  Discoveries,  Inventions,  Patent  Disclosures 


No  new  patent  disclosures. 


9.0  Honors  and  Awards 


None. 


The  Optimal  Design  of  Structures  Incorporating  ‘Smart  Materials’ 


9  of  9 


Optimal  distribution  of  material  properties  for 
an  elastic  continuum  with  structure-dependent 

body  force. 

Sergio  Turteltaub  and  Peter  Washabaugh 


Department  of  Aerospace  Engineering 
The  University  of  Michigan 
Ann  Arbor,  MI  48109 


Abstract 

The  simultaneous  optimization  of  material  properties  and  struc¬ 
tural  layout  for  an  elastic  continuum  is  formulated  and  analyzed.  The 
objective  is  to  obtain  the  maximum  structural  stiffness  for  prescribed 
surface  loads  and  displacements,  taking  into  account  a  body  force 
that  depends  on  the  structural  layout.  Optimization  with  an  account 
of  self-weight  or  centrifugal  forces  are  examples  of  these  type  of  prob¬ 
lems.  Arbitrary  elasticity  tensor  fields  are  considered  as  the  problem’s 
variable  and  necessary  conditions  satisfied  by  the  solution  are  estab¬ 
lished.  The  use  of  a  spectral  decomposition  of  the  elasticity  tensor  is 
emphasized  since  it  provides  a  simple  geometrical  interpretation.  Typ¬ 
ical  examples  which  illustrate  the  effects  of  the  structure-dependent 
body  force  are  analyzed.  It  is  found  that  a  commonly  used  isoperi- 
metric  restriction  (known  as  the  resource  constraint)  is  not  necessar¬ 
ily  active  and  that  the  optimal  structure  is  locally  stiffer  in  areas 
where  the  body  force  and  the  displacement  field  have  opposite  di¬ 
rections  and  locally  weaker  otherwise.  This,  interestingly,  leads  to  a 
non-symmetrical  distribution  of  material  properties  even  for  prismatic 
bodies  under  anti-symmetric  surface  loads. 
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1  Introduction. 

The  problem  of  characterizing  the  maximum  structural  stiffness  of  a  linearly 
elastic  continuum  structure  has  been  extensively  analyzed  using  different 
techniques.  It  has  been  recognized  that  it  is  possible  to  find  simultaneously 
the  optimal  material  properties  and  structural  layout.  One  way  to  achieve 
this  is  simply  by  enlarging  the  space  of  variables  in  order  to  include  all 
possible  structures,  i.e.,  by  considering  completely  general  elasticity  tensor 
fields.  This  approach  is  known  as  the  free-material  optimization  method  and 
it  relies  on  the  use  of  elasticity  tensor  fields  which,  a  priori,  correspond  to 
anisotropic  and  inhomogeneous  materials  (see  Bendsoe  et  al.  (1994)).  Other 
methods  consider  a  special  class  of  anisotropic  materials,  namely  composites 
obtained  by  combining  two  distinct  homogeneous  materials,  although  the 
problem  requires  relaxation  in  order  to  include  optimal  composites.  In  that 
case,  usually  assuming  the  existence  of  some  type  of  periodic  microgeometry, 
homogenization  theory  is  used  to  compute  the  effective  properties  at  a  macro¬ 
scopic  level.  Parameters  related  to  the  microgeometry  are  used  as  variables 
(see,  e.g.,  Bendsoe  (1995)  for  a  comprehensive  review  of  both  methods).  The 
advantage  of  the  free-material  optimization  method  is  that  issues  related  to 
homogenization  do  not  need  to  be  taken  into  account.  Once  the  optimal  ma¬ 
terial  properties  have  been  identified  (at  a  macroscopic  level)  one  can  solve,  if 
desired,  an  inverse  problem  in  order  to  obtain  at  each  point  a  microgeometry 
that  provides  the  best  match  to  the  optimal  properties  as  shown  by  Sigmund 
(1994)  (see  also  Milton  and  Cherkaev  (1995)).  The  purpose  of  the  present 
analysis  is  to  include,  within  the  free-material  formulation,  a  body  force  that 
depends  on  the  optimization  variable  (i.e.,  on  the  elasticity  tensor).  The 
motivation  for  this  extension  is  to  analyze  cases  where  a  non-negligible  in¬ 
ertial  force  acts  on  the  structure.  Typical  examples  are  structural  elements 
rotating  at  a  relatively  large  angular  speed  or  large  structures  whose  own 
weight  becomes  a  relevant  factor  in  the  analysis.  Early  work  on  problems 
with  a  structure-dependent  body  force  include  applications  where  the  ob¬ 
jective  is  to  distribute  a  given  amount  of  material  in  an  elastic  structure  in 
order  to  match  a  desired  natural  frequency  (see,  e.g.,  Prager  (1974)  and  the 
references  therein).  A  more  recent  review  on  the  subject  was  given  by  Olhoff 
(1987).  Here,  however,  the  objective  is  to  maximize  the  structural  stiffness 
taking  into  account  a  variable  body  force  but  no  restrictions  are  placed  upon 
natural  frequencies. 

The  use  of  a  formulation  based  on  a  spectral  decomposition  of  the  elas- 
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ticity  tensor  is  emphasized  since  it  provides  a  simple  geometrical  interpre¬ 
tation  of  the  results  in  a  fourth-order  tensor  space.  It  will  be  shown  that 
the  optimized  material  has  two  main  components:  a  term  which  provides 
the  optimal  stiffness  and  a  term  which  provides  the  required  stability  against 
possible  perturbations  of  the  prescribed  loads.  In  this  analysis,  it  is  assumed 
that  the  magnitude  of  the  elastic  moduli  and  the  material’s  mass  density  are 
correlated.  Even  though  in  general  there  is  no  a  priori  relation  between  these 
quantities,  some  models  used  for  cellular  solids  exhibit  a  functional  relation 
between  the  elastic  moduli  (see,  e.g.,  Gibson  and  Ashby  (1997)).  From  a  dif¬ 
ferent  point  of  view,  a  relation  between  mass  density  and  elastic  properties 
can  be  assumed  if  one  has  in  mind  using,  a  posteriori,  a  method  like  the  one 
proposed  by  Sigmund  (1994).  In  that  case,  a  microgeometry  composed  of 
weak  (essentially  void)  and  strong  materials  is  used  to  match  given  effective 
properties.  Bearing  in  mind  that  the  strong  material  is  fixed  (and  so  is  its 
mass  density),  the  resulting  volume  fraction  of  the  strong  material  is  directly 
related  to  the  mass  density  of  the  composite  material.  Thus,  in  this  restricted 
sense,  one  can  establish  a  relation  between  elastic  moduli  and  mass  density 
and  it  is  assumed  that  this  relation  is  strictly  monotonic. 

The  outline  of  the  paper  is  as  follows:  in  Section  2,  the  optimization  prob¬ 
lem  is  formulated  and  under  some  specific  assumptions  a  simplified  version 
of  it  is  derived  by  analyzing  the  restrictions  imposed  at  a  local  level  (for  a 
continuum)  by  the  global  structural  requirements.  Optimality  conditions  for 
the  problem  are  developed  in  Section  3  and  specific  examples  with  an  inertial 
body  force  are  treated  in  Section  4.  A  discussion  and  concluding  remarks 
follow  in  the  last  section. 

2  Formulation  of  the  problem. 

2.1  Notation. 

As  a  general  scheme  of  notation,  scalar  quantities  are  represented  by  italicized 
normal-face  letters,  vectors  and  points  in  a  three-dimensional  Euclidean  space 
by  bold-face  lower  case  letters  (except  for  the  stress  tensor  cr  and  strain 
tensor  e),  second  order  tensors  by  bold- face  upper  case  letters  and  fourth 
order  tensors  by  bold-face  upper  case  italicized  letters.  Throughout  this 
analysis,  coordinate-free  notation  is  employed.  To  ease  the  transition  to 
indicial  notation,  the  following  definitions  are  referred  to  a  three-dimensional 
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Cartesian  basis  (indices  range  in  {1,2,3},  6ij  is  Kronecker’s  delta  and  the 
summation  convention  is  used): 


a  •  b  =  g,-6,'  , 
(Ab);  =  Aijbj  , 


A  •  B  =  AijBij  , 

II A||  =  VaTa  , 

(AB)ij  =  AikBkj  , 
(a  0  b),-j  =  Uibj  , 

(I h  =  . 


CD  —  CijkiDijki  , 

\\C\\  =  y/C~C  , 

{CA)ij  =  CijkiAki  , 

(A  0  B)iyjfc/  =  AijBu  , 
(I)ijki  =  $(8ik5ji  +  SiiSjk)  . 


In  general,  when  the  meaning  is  clear  by  the  context,  a  scalar,  vector  or  tensor 
field  defined  in  a  domain  Cl  and  its  value  at  a  point  x  will  be  denoted  by  the 
same  letter,  e.g.,  A  =  A(x).  Similarly,  a  scalar,  vector  or  tensor  function 
of  a  scalar,  vector  or  tensor  field  will  be  denoted  by  the  same  symbols,  e.g., 
ip  =  <p{A)  =  <p(A(x)),  when  it  is  clear  by  the  context  that  <p  is  not,'  for 
example,  a  functional. 


2.2  Preliminaries. 

Consider  a  linearly  elastic  body  that  occupies  a  given  domain  Cl.  Suppose 

that  the  prescribed  tractions  t  and  displacements  u  are 

» 

t(x)=t(x)  xedClt  \  m 

u(x)  =  0  x  6  dClu  j  K  1 

where  3n£u5flu  =  dfl  is  the  boundary  of  Cl.  Only  mixed  boundary  conditions 
of  the  form  (1)  are  considered  here.  Let  e  be  the  strain  tensor  field  which, 
viewed  as  a  function  of  u(x),  is  given  by 

e  =  e(u)  =  e(u(x))  =  ^  (Vu(x)  +  Vu(x)T)  ,  Vx  €  Cl  . 

Let  C  be  a  (fourth-order)  elasticity  tensor  field.  The  stress  tensor  field  a  is 
related  to  the  strain  tensor  via  the  constitutive  law 


cr(x)  =  C(x)s(x)  . 

As  a  matter  of  terminology,  the  word  “structure”  should  be  understood  in 
the  present  analysis  to  refer  to  a  distribution  of  material  properties  over  the 
given  domain  Cl,  i.e.,  to  the  field  C(x)  ,  x  €  Cl.  By  “optimal  distribution”  it 
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is  meant  that  a  specific  structure  (i.e.,  distribution  of  C(x))  is  sought  based 
on  maximization  of  the  overall  stiffness  of  the  structure.  A  formal  statement 
of  this  problem  is  given  in  Section  2.4. 

To  obtain  the  optimal  distribution  of  material  properties,  the  basic  idea 
is  to  let  the  elasticity  tensor  field  be  variable.  It  is  worth  noting  that  one 
can  accomplish  two  seemingly  distinct  objectives  with  this  approach.  Firstly, 
the  minimizing  field  C o  provides  information  on  the  mechanical  properties 
of  a  material  at  each  point  (i.e.,  the  material  optimization).  Secondly,  since 
points  where  the  material  properties  are  “small”  are  interpreted  as  holes  (in 
a  sense  to  be  defined  later),  then  C0  also  provides  the  shape  of  the  struc¬ 
ture  within  the  domain  Q  by  specifying  the  location  of  holes.  Hence,  the 
shape  optimization  (structural  layout)  is  obtained  as  a  by-product  of  the 
material  optimization.  This  ability  to  perform  simultaneously  optimization 
of  material  properties  and  structural  layout  is  one  of  the  main  benefits  of  the 
free-material  optimization  method  and  its  simplicity  (compared  to  homoge¬ 
nization  techniques)  provides  an  additional  advantage. 

2.3  Body  force  and  space  of  admissible  variables. 

The  space  of  admissible  optimization  variables  is,  in  this  case,  a  subspace  of 
fourth-order  tensor  fields  defined  in  fl.  For  practical  reasons,  both  physical 
and  mathematical,  some  additional  restrictions  on  the  space  of  optimization 
variables  are  required.  To  introduce  these  restrictions  it  is  convenient  to 
use  a  spectral  decomposition  of  the  elasticity  tensor.  Since  the  admissible 
elasticity  tensors  are  symmetric  (major  symmetry),  they  can  be  represented 
at  each  point  x  in  terms  of  their  real  eigenvalues  (p  =  1,. . .  ,6)  and  unit 
eigentensors  (2-tensors),  i.e.,  the  spectral  decomposition  is  given  by 

6 

c  =  ^2  °V >  (2) 

M=1 

where  >  0  and  {AM}®=1  forms  an  orthonormal  basis  for  the  space  of 
symmetric  2-tensors.  Recall  that  for  a  general  (anisotropic)  material,  the 
eigentensors  contain  information  about  the  material  behavior  (via  “angles” 
in  a  tensor  space)  but  the  “true”  stiffness  information  is  provided  by  the  six 
principal  values  of  the  elastic  moduli  (i.e.,  eigenvalues  a,,).  For  example,  for 
an  isotropic  material  there  are  two  distinct  eigenvalues:  3/c  (multiplicity  one) 
and  2p  (multiplicity  five),  where  k  is  the  bulk  modulus  and  p  is  the  shear 
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modulus  (see,  e.g.,  Knowles  (1995)).  It  is  also  important  to  note  that  the 
result  of  this  decomposition  is  that  the  design  variables  are  now  explicitly 
the  eigenvalues  and  eigenvectors. 

The  fact  that  the  eigenvalues  are  strictly  positive  means  that  C  is  posi¬ 
tive  definite.  To  satisfy  this  requirement,  it  is  assumed  that  all  eigenvalues 
are  bounded  below  by  the  same  constant  am.  Additionally,  an  upper  bound 
ccm  is  imposed  on  each  (same  constant  ccm  for  all  ap).  Since  {o:^}®=i 
are  the  principal  stiffnesses,  the  upper  bound  is  seen  as  a  limitation  of  real 
materials:  the  elastic  moduli  cannot  exceed  the  values  of  the  stiffest  material 
known  in  nature  (or,  in  fact,  of  a  predetermined  material  that  would  be  used 
to  reinforce  the  structure).  Finally,  a  global  upper  bound  on  the  eigenvalues 
(known  as  a  resource  constraint )  is  also  imposed.  To  this  end,  consider  a 
positive-valued  scalar  function  (p  of  C  which  is  referred  to  as  the  resource 
constraint  density.  In  the  work  of  Bendspe  et  al.  (1994),  a  resource  constraint 
density  equal  to  the  trace  of  C  was  used,  i.e.,  < p(C)  =  Following  a 

similar  approach,  the  resource  constraint  density  considered  here  is  assumed 
to  be  independent  of  the  orientation  of  the  principal  directions  of  the  elastic¬ 
ity  tensor  and  dependent  on  the  eigenvalues  of  C  via  their  maximum  only, 
i.e., 


0(C)  =  a,  (3) 

where  a  is  the  maximum  principal  stiffness,  i.e., 

a  =  q(x)  =  maoc  (ap(x)}  (4) 

l<ji<6 

Observe  that  (3)  is  not  the  only  choice  since  (p  could  depend  on  (the 
orientation  of  material  properties)  which  could  be  interpreted,  e.g.,  as  a  con¬ 
straint  imposed  by  manufacturing  requirements.  In  terms  of  a,  the  resource 
constraint  is  expressed  by  the  following  isoperimetric  inequality: 


/  a(x)du  <  R  , 

J  n 

where  the  left  hand  side  is  the  total  resource  and  the  given  positive  con¬ 
stant  R  represents  an  upper  bound  on  the  resource.  The  main  role  of  the 
resource  constraint  is  to  rule  out  some  trivial  solutions.  A  typical  trivial  case 
occurs  when  the  body  force  b  is  zero  and  the  isoperimetric  inequality  is  not 
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enforced:  the  stiffest  structure  is  obtained  when  the  optimization  variable  co¬ 
incides  with  the  local  upper  bound,  i.e.,  it  coincides  with  the  stiffest  material 
properties  everywhere. 

For  future  use,  all  the  requirements  on  the  space  of  admissible  optimiza¬ 
tion  variables  can  be  collected  via  the  following  sets: 

*Sap  =  j{a„}J=si  |  ccm  <  aM(x)  <  ocM  ,  J^adv  <R,  Vx  e  Q 

and 

<$a„  =  I-^-p  '  =  ^ ~  Aft  >  ®  ~  1 

where  a  is  defined  by  (4)  and  6^u  is  Kronecker’s  delta.  The  symmetry  of 
each  Af,  reflects  the  minor  symmetries  of  C.  Now,  it  is  assumed  that  there 
exist  functional  relations  between  the  principal  elastic  moduli  (eigenvalues 
dp)  and  the  mass  density  p,  i.e.,  =  aM(p).  To  motivate  this  assumption 

from  a  physical  point  of  view,  one  could  refer,  for  example,  to  models  used 
for  cellular  solids  (see  Gibson  and  Ashby  (1997)).  In  that  context,  it  has 
been  found  experimentally  that,  in  the  linearly  elastic  range,  the  material 
properties  scale  with  some  power  of  the  mass  density  (e.g.,  the  shear  modulus 
for  elastomeric  foams  is  p  ~  p2,  etc.).  Similar  relations  are  sometimes  used 
in  some  simple  models  for  composite  materials  yvhere  the  material  properties 
depend  on  the  volume  fraction  of  each  component  (hence  they  depend  on  the 
relative  mass  densities).  Furthermore,  assuming  that  these  relations  can  be 
inverted,  then  the  mass  density  could  be  expressed  as  a  function  of  one  of  the 
eigenvalues  of  C;  in  particular,  p  =  p(a),  where  a  is  the  maximum  principal 
stiffness.  Thus,  assuming  the  existence  of  such  a  model,  an  (inertial)  body 
force  can  be  viewed  as  a  function  of  a,  i.e., 

b  =  b(p(a))  =  b(a)  .  (5) 

With  this  interpretation,  the  resource  constraint  can  also  be  seen  as  a  restric¬ 
tion  on  the  total  mass  of  the  structure.  However,  it  is  important  to  mention 
that  the  present  analysis  is  limited  to  cases  where  a  relation  such  as  (5)  can 
be  identified. 

Let  W  be  the  work  done  by  the  external  forces,  i.e., 

W[u,b(a)]=  [  b(a(x))  •  u(x)dv  +  J  t(x)  •  u(x)da  . 

Jn  h  n, 


(6) 


Optimal  distribution  of  material  properties.  (Rev.  form) 


8 


Let  V  be  a  suitably  chosen  space  of  kinematically  admissible  displacement 
fields.  For  given  boundary  conditions  (1)  and  for  a  specific  (i.e.,  “frozen”) 
elasticity  tensor  field  C,  the  displacement  field  u  that  satisfies  the  equilibrium 
equation  is  the  (unique)  element  of  the  set  Sc  defined  as  follows: 

Sc  =  ju  6  V  |  J  e(u)  •  Ce(v)du  =  W(v,  b(a)]  Vv  e  V  j  .  (7) 

The  work  W  given  by  (6)  is  viewed  as  a  global  measure  of  the  stiffness  of 
a  structure  and  takes  into  account  the  prescribed  domain  Q  and  boundary 
conditions.  The  stiffest  structure  corresponds  to  the  one  which  minimizes  W 
for  fixed  boundary  conditions  (1)  when  the  elasticity  tensor  field  C  is  taken 
as  variable. 


2.4  Optimization  problem:  maximum  structural  stiff¬ 
ness. 

The  optimization  problem  is  stated  as:  for  a  linearly  elastic  material  occu¬ 
pying  a  given  domain  Q  and  for  boundary  conditions  given  by  (1),  find  the 
minimizer  C  =  £®=i  ot^  A„  ®  Ap  of  the  following  expression: 


min  min  HTu,  b(a)l 
{Am}€«Sap 


(pi) 


where  u  €  Sc  (hence  the  admissible  field  u  is  the  equilibrium  solution  for  a 
given  field  C)  and  at  is  given  by  (4).  The  two  minimization  parts  reflect  a 
decomposition  of  the  optimization  problem  in  terms  of  two  sets  of  variables 
(eigenvalues  {a*,}  and  eigentensors  {A^}  of  C).  The  sets  <Sa>1  and  <SA„ 
correspond  to  the  constraints  on  the  optimization  variables.  For  the  solution 
u  of  the  elastostatic  problem  (i.e.,  u  e  Sc),  the  dependence  of  W  on  C ,  as 
given  by  (6),  is  two-fold:  u  depends  implicitly  on  the  constitutive  law  and,  by 
assumption,  the  body  force  depends  explicitly  on  the  maximum  eigenvalue  of 
C.  Formally,  W  also  depends  on  fi  and  the  prescribed  boundary  conditions, 
but  these  are  considered  fixed. 

The  problem  (PI)  can  be  simplified  considerably  by  virtue  of  a  saddle 
point  theorem  as  proved  in  Bendsoe  et  al.  (1994)  (see  also  Jog  et  al.  (1993) 
and  Lipton  (1994)).  This  is  an  established  result  in  the  optimization  of 
structures.  Nonetheless,  in  order  to  illustrate  the  explicit  role  played  by 
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the  eigenvalues  and  eigentensors  of  C  introduced  in  (2)  and  to  show  that 
the  saddle  point  theorem  remains  unaffected  by  a  structure-dependent  body 
force  of  the  form  (5),  it  is  useful  to  revisit  this  result.  More  significantly, 
the  analysis  presented  here  provides  an  interpretation  of  the  solution  of  the 
optimization  problem  in  terms  of  a  positive  definite  material,  as  opposed 
to  previous  solutions  which  were  given  in  terms  of  semi-definite  materials. 
To  this  end,  the  admissible  displacement  fields  in  (Pi)  are  characterized 
as  minimizers  of  the  potential  energy  (instead  of  enforcing  the  equilibrium 
equations).  The  potential  energy  is  given  by  II  =  U  —  W,  where  U  = 
5  /n  e '  Cedu  is  the  total  strain  energy.  The  minimum  value  of  the  potential 
energy  is  II.  =  —  \W,,  hence,  reversing  the  sign  of  the  objective  functional, 
the  problem  (Pi)  can  be  expressed  alternatively  as:  find  the  maximizing 
eigenvalue  fields,  eigentensor  fields  and  minimizing  displacement  field  of  the 
following  expression: 

max  max  minII[aW)  A„;  v]  ,  (P2) 

{a„}es0„  {AM}esA„  v6V 

where,  since  e  •  Ce  =  aM(e  •  A M)2,  the  potential  energy  is 


The  minimization  part  in  (P2)  corresponds  to  the  elastostatic  problem  for  a 
given  elasticity  tensor  field.  In  view  of  the  saddle  point  theorem  mentioned 
above,  the  two  inner  problems  in  (P2)  can  be  interchanged  (the  theorem  can 
be  applied  because  it  is  assumed  that  b  is  not  a  function  of  Ap).  It  is  worth 
noting  that  the  outermost  maximization  problem  cannot  be  interchanged 
with  the  innermost  minimization  problem  since  they  provide  different  solu¬ 
tions.  After  interchanging  the  two  inner  problems  in  (P2),  the  innermost 
problem  becomes,  for  given  v  and  {ap}®=1,  a  local  algebraic  problem,  i.e., 
find  the  maximizers,  at  each  point  x  €  12,  in  the  following  expression: 

max  zMe(v) '  AJ2  ■  (9) 

{A„}€SA(,  " 
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To  provide  a  “geometrical”  interpretation  of  the  problem  notice  that,  since 
e  ■  Ce  =  (e  ®  e)  •  C,  then  (9)  is  equivalent  to 

max 

where  the  maximization  is  carried  out  for  six  4-tensors  of  the  form  A,,  ®  Ap 
and  the  problem  is  linear  since  the  goal  is  to  maximize  the  scalar  product  of 
“vectors”  in  a  4-tensor  space,  where  both  and  e®  e  are  fixed.  The 

objective  function  in  the  inner  problem  (10)  (or  (9))  admits  the  following 
trivial  bound: 

6 

|(e®e)  -^MAm®  K)  <  5ma^6{orp}||e®e||  =  ^ae-e  .  (11) 

n=\ 

This  upper  bound  is  independent  of  {Ap}  and  depends  on  the  local  values 
of  {aM}  only  via  their  maximum  a.  To  reach  this  bound,  one  could  choose 
the  eigentensor  corresponding  to  the  maximum  eigenvalue  as 

A  =  s  =  — 

Nl’ 

while  the  other  eigentensors  can  be  chosen  arbitrarily  as  long  as  they  form  an 
orthonormal  basis.  This  choice  can  be  made  regardless  of  the  multiplicity  of 
the  maximum  eigenvalue  (i.e.,  one  of  the  eigentensors  corresponding  to  a  can 
be  defined  as  above).  With  this  specific  choice  the  upper  bound  is  reached. , 
thus  providing  an  optimal  solution  for  the  set  {A^,},  valid  for  any  point 
x  €  Q,  in  terms  of  the  (kinematically  admissible)  displacement  field  v.  An 
explicit  expression  for  the  eigentensors  different  than  k  is  not  required  for  the 
present  purposes.  From  (8),  (11)  and  the  definition  of  5Q(J,  one  can  see  that 
if  the  upper  bound  is  reached,  then  the  problem  (P2)  depends  only  on  the 
maximum  eigenvalue  a.  In  that  case,  the  eigenvalues  smaller  than  a  can  be 
chosen  arbitrarily  as  long  as  they  belong  to  SQtl.  This  problem  has  no  unique 
solution,  however  a  convenient  choice  is  to  set  these  eigenvalues  to  the  lower 
bound  am  in  order  to  guarantee  that  a  is  the  maximum.  If  the  multiplicity 
of  a  is  greater  than  one,  then  the  above  choice  corresponds  to  reducing  its 
multiplicity  to  one.  The  essential  point  is  that  any  admissible  choice  provides 
the  same  final  result  (i.e.,  upper  bound).  The  optimal  material  can  now  be 


|(e(v)  ®  e(v))  ®  Ap) 


(10) 
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expressed  as 

C0  =  C0(a,v)  =  (a  -  -  +  amI  .  (12) 

The  coupling  between  the  optimization  and  elastostatic  problems  is  reflected 
in  the  fact  that  the  elasticity  tensor  depends  on  the  strain  field.  It  is  noted  in 
passing  that  Co  is  an  orthotropic  material  with  symmetry  directions  aligned 
(locally)  with  the  principal  directions  of  the  strain  tensor:  to  see  this,  sup¬ 
pose  that  Q  is  a  proper  orthogonal  second  order  transformation  that  com¬ 
mutes  with  e,  hence  Q  and  e  have  the  same  principal  directions  in  an  Eu¬ 
clidean  space.  Since  Co  transforms  as  x  or/1QA^QT®  QApQr  and  since 
QeQT  =  e  then  it  follows  from  (12)  that  Q  belongs  to  the  material  sym¬ 
metry  group  of  Co-  The  orthotropy  of  the  optimal  material  is,  in  fact, 
a  well-known  result  (see,  e.g.,  Pedersen  (1989),  Cowin  (1994)).  However, 
the  specific  form  (12)  has  a  new  ingredient  compared  to  the  one  proposed 
by  Bendsoe  et  al.  (1994),  i.e.,  the  term  omJ,  which  provides  the  required  sta¬ 
bility  of  the  material.  .Recently,  Taylor  (1997)  proposed  a  general  framework 
where,  among  other  things,  terms  such  as  amI  can  be  easily  introduced. 
The  essential  term  in  the  material  (12)  is  a  4-tensor  perpendicular  projec¬ 
tion  onto  the  2-tensor  space  spanned  by  the  strain  tensor  (which  can  be 
interpreted  as  an  optimal  “reinforcement,”  characterized  by  the  maximum 
principal  stiffness  a).  For  illustration  purposes,  the  material  properties  of 
the  orthotropic  material  can  also  be  expressed  in  terms  of  the  (three)  Young 
moduli  E{,  the  (three  independent)  Poisson’s  ratios  and  (three)  shear 
moduli  Gij.  Let  i  =  1,2,3,  be  the  principal  values  of  the  (unit)  strain 
tensor  i:  =  £/||e||;  for  an  arbitrary  strain  tensor  7,  the  corresponding  stress 
is,  from  (12),  a  =  C07  =  (a  —  am)(£  ~l)e  +  Qm7-  Hence,  the  elastic  moduli 
referred  to  a  Cartesian  basis  aligned  locally  with  the  principal  directions  of 
£(v)  are  given  by 


e?)a  1 


(a  -  Qm)g«gj 
am  +  (1  -  £?)a  ’ 


Gij  —  2^r' 


Moreover,  the  stress  field  associated  with  the  optimal  strain  field  £(v)  takes 
the  following  simple  form: 


<x  =  C0e  =  <ye  . 


(13) 


By  optimal  strain  field  it  is  meant  that  £  is  the  strain  field  in  the  elastostatic 
problem  after  the  material  C  has  been  optimized  analytically  with  respect  to 
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the  set  {Ap}.  Although  the  elasticity  tensor  given  by  (12)  is  anisotropic,  its 
restriction  to  the  optimal  strain  field,  given  by  (13),  corresponds  formally  to 
an  isotropic  material  with  material  coefficients  3k  =  2p  =  a  or,  equivalently, 
to  v  —  0  and  E  =  a  where  v  is  Poisson’s  ratio  and  E  is  Young’s  modulus. 
This  formal  interpretation  is  valid  since  the  optimization  and  elastostatic 
problems  are  being  solved  simultaneously,  hence  the  optimal  material  prop¬ 
erties  are  coupled  to  the  optimal  displacement  field  (this  is  referred  to  as  a 
self-adaptive  material,  although  it  is  a  mathematical  rather  than  a  physical 
property). 

In  view  of  (12),  the  problem  (Pi)  can  be  simplified  in  terms  of  the  scalar 
field  a  (maximum  principal  stiffness),  i.e.,  find  the  minimizer  a  of 

min  W[u,  b(a)]  ,  (P3) 

o€Sa 

where  u  €  £a, 

Sa  =  |u  €  V  |  /  ae(u)  •  e(v)du  =  W(v,  b(o:)]  Vv  €  V  j  (14) 

and 

<SQ  =  ja|  om  <  q(x)  <  Qa  {  ,  J  adu  <  R  ,  Vx  £  ft  j 

At  points  where  the  maximum  principal  stiffness  a  is  equal  to  its  lower  bound 
Qm,  the  interpretation  is  that  no  reinforcement  is  required.  The  formulation 
(P3)  provides  a  significant  simplification  from  (Pi)  since  the  minimization  is 
carried  out  for  a  single  scalar  field  as  opposed  to  a  tensor  field.  For  given  a, 
the  local  form  of  the  elastostatic  problem  for  the  optimal  strain  e  at  points 
where  a  and  e  are  smooth  is,  from  (14), 

div(a£r)  +  b(a)  =  0  in  f2  ,  'j 

crn  =  aen  =  t  on  dQt  ,  ?  (15) 

u  =  0  on  dflu  ,  J 

where  n  is  the  outward  unit  normal  vector  to  the  boundary  dfl.  Observe 
that  the  for  shape  optimization,  “holes”  are  identified  via  a  limit  process 
when  Om/ctM  — »  0,  for  which  uniqueness  is  preserved.  However,  in  the  limit, 
the  material  becomes  semi-definite,  which  is  in  fact  one  of  the  limitations 
of  the  “true”  single-loading  shape  optimization  problem.  Hence,  after  the 
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solution  to  (P3)  has  been  identified  (optimal  structure)  and  is  considered  to 
be  fixed,  am  needs  to  be  bounded  away  from  zero  by  an  amount  comparable 
to  possible  perturbations  on  the  surface  loads  in  order  to  have  a  “robust” 
structure. 


3  Optimality  conditions. 


The  optimality  conditions  correspond  to  a  set  of  necessary  relations  satisfied 
locally  by  the  minimizer  of  (P3)  and  can  be  used  for  numerical  methods  or 
to  solve  some  simple  problems  in  closed  form  as  will  be  shown  in  Section  4. 
To  obtain  these  conditions,  consider  the  Lagrangian  L  given  by 


L[or,  Am,  \m,  A]  =  W0[u,  b(a)]  +  /  [Am(o;m  -  a)  -  Xm(q  -  qw)]  dv 

Jo 

,  +A{/„ad“-d . 


where  Am  =  Am(x),Aw-  =  Aw(x)  are  Lagrange  multipliers  associated  with 
the  upper  and  lower  bound  constraints  respectively  and  A  is  the  (constant) 
multiplier  corresponding  to  the  resource  constraint.  As  mentioned  before, 
the  displacement  field  u  is  viewed  as  a  function  of  a  since  it  can  be  obtained 
from  (14).  Hence,  in  terms  of  the  Lagrangian  L ,  the  problem  (P3)  can  be 
expressed  as:  find  the  maximizers  Am,  \m,  A  and  minimizer  a  in  the  following 
expression: 


max  minLfa;  Am,  A^,  A]  . 

Am>0,  am>o,  a>o  a 

The  gradient  of  L  with  respect  to  a  can  be  computed  as  follows:  consider  a 
variation  Sec  which  induces  variations  <5u  and  <5b.  Let  V  be  the  gradient  of 
L,  therefore,  from  (5),  (6),  (12)  and  (16),  it  follows  that 

6L[5oc\  =  f  L'Sadv  =  j  (<5b  •  u  +  b  •  6u )  dv  +  f  t  •  <5uda 
J  n  J  n  Jdot 

-  /  (Am  -  Am  -  A)  Socdv  . 

Jo 


(17) 
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To  complete  the  calculation,  one  can  take  variations  in  (14),  then  integrate 
by  parts,  choose  v  =  u,  and  use  (15)  and  (17)  to  get 

f  L'5adv  =  f  25b  •  udu  —  f  5a  (e(u)  •  e(u)  +  Am  -  Xm  —  A)  dv  . 

Jn  J  n  J  n  ' 

Since  5b  =  b'5a,  where  b'  is  the  gradient  of  b  with  respect  to  a,  then, 
assuming  enough  differentiability,  the  local  form  of  the  gradient  is 

U  =  2b'  •  u  -  e(u)  •  e(u)  -  Am  +  +  A  .  (18) 

Therefore,  the  optimality  conditions  (Karush-Kuhn-Tucker  conditions)  are, 
Vx  6  0, 

U  =  0  , 

Am(<am  -  a)  =  0  ,  Am  >  0  , 

—  aM)  =  0  i  ^  0  . 

A  {/„  adv  -  J?}  =  0  ,  A  >  0  . 

To  interpret  the  optimality  conditions  one  can  use  the  following  domains: 
nm  =  (x  €  Q  |  or(x)  =  Qm}  , 

=  {x  6  ft  |  q(x)  =  aM }  ,  (20) 

fii  =  (x  e  n  I  Qm  <  q(x)  <  aM)  . 

Hence,  since  Am  =  0  for  x  6  Qm,  =  0  for  x  €  Qm  and  Am  =  Xm  =  0  for 
x  €  Qi,  the  optimality  conditions  (19)  become 

A  =  e(u)  •  e(u)  —  2b'  •  u  , 

A  —  Am  =  e(n)  ■  e(u)  -  2b'  •  u  ,  Am  >  0  , 

A  +  A m  —  s(u)  •  e(u)  -  2b'  •  u  ,  A^-  >  0  , 

A  {/n  adv  —  R}  =  0  ,  A  >  0  , 

If  b  =  b'  =  0  then  one  can  prove  that  A  >  0  (and  hence  that  the  resource 
constraint  is  active  (i.e.,  satisfied  as  an  equality;  see  Bendsese  et  al.  (1994)). 
Clearly,  to  have  a  non-trivial  solution  in  the  case  when  the  resource  constraint 
is  active,  the  upper  bound  R  in  the  resource  constraint  has  to  satisfy  am|Q|  < 
R  <  a/vfini-  However,  if  b'  ^  0,  it  is  possible  that  A  =  0  and  the  problem 
(P3)  has  a  non-trivial  solution.  This  case  will  be  illustrated  by  an  example  in 


x€  , 

X  6  om  ’  >  (21) 

x  €  Qm  , 
x  €  Q  . 


Optimal  distribution  of  material  properties.  (Rev.  form) 


15 


Section  4.  Furthermore,  recall  that  it  was  assumed  that  b  is  a  monotonically 
increasing  function  of  oc  hence,  in  view  of  (18),  if  b  •  u  >  0,  the  local  effect  of 
the  term  2b'  •  u  is  to  increase  the  value  of  the  gradient  of  the  Lagrangian  L. 
Conversely,  if  b  •  u  <  0  (hence  b'  •  u  <  0),  the  local  effect  of  the  term  2b'  •  u 
is  to  decrease  V .  Therefore,  compared  to  the  case  when  the  body  force  zero, 
the  optimal  structure  is  locally  weaker  (smaller  maximum  principal  stiffness 
a)  if  b  •  u  >  0  and  locally  stiffer  (greater  a)  otherwise.  This  effect  will  be 
illustrated  by  a  three-dimensional  example  in  next  section. 


4  Example:  the  optimal  rotating  structure. 

As  an  example  of  a  structure-dependent  body  force,  consider  the  problem  of 
finding  the  stiffest  structure  occupying  a  domain  fi  and  rotating  at  a  constant 
angular  velocity  u.  Let  {er,  eg,  ez }  be  a  cylindrical  orthonormal  basis  where 
ez  is  aligned  with  the  axis  of  rotation  and  eT,e$  are  fixed  with  respect  to 
fi  (er  points  in  the  outward  radial  direction).  For  illustration  purposes, 
suppose  that  the  optimization  is  carried  out  with  a  material  for  which  the 
mass  density  p  is  related  to  the  maximum  principal  stiffness  a  linearly,  i.e., 

p  —  KQr  , 

where  k  is  a  positive  constant.  It  is  worth  noting  that  for  nonlinear  functions 
p  the  corresponding  problem  is  nonlinear  and  might  not  have  a  solution.  The 
inertial  force  can  be  modeled  in  the  elastostatic  problem  as  the  body  force 

b(a)  =  u?rpeT  =  u;2r/caer  ,  (22) 

where  r  is  the  orthogonal  distance  from  a  given  point  in  fi  to  the  axis  of 
rotation.  To  obtain  some  insight  in  the  problem,  it  is  possible  to  solve  ana¬ 
lytically  a  one-dimensional  example.  To  this  end,  consider  a  prismatic  region 
fi  of  length  /  and  cross-sectional  area  a.  Suppose  that  the  displacement  field 
and  the  maximum  principal  stiffness  a  depend  on  r  only,  r  €  [0,/].  As  an 
ansatz,  assume  that  the  domains  defined  by  (20)  are  such  that  fi m  =  (0,  ri], 
fi,-  =  [rl,r2],  fim  =  [r2 ,  Z]  and  that  the  resource  constraint  is  not  active 
(i.e.,  A  =  0).  It  will  be  shown  a  posteriori  that  these  assumptions  are  ap¬ 
propriate  to  characterize  solutions  with  boundary  conditions  u(0)  =  0  and 
ct  =  a(l)u'(l )  =  oo  >  0,  where  u(r)  is  the  radial  displacement.  The  optimality 
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conditions  (21)  become 

A,w(r)  =  ( u'(r ))2  —  2 ku2tu(t)  A iVf(r)  >0  0  <  r  <  77 

(u'(r))2  —  2 kuttu(t)  =  0  77  <  T  <  r2  (23) 

Am(r)  =  -(u'(r))2  +  2KU2ru(r)  A m(r)  >0  r7  <  r  <  l 

and  the  balance  of  linear  momentum  is,  from  (15)  i , 

(ool1)'  +  nu7ra  =0  0  <  r  <  l  .  (24) 

Observe  that,  since  a  =  cxm  >  0  in  ft.v  and  a  =  am  >  0  in  ftm,  then  (24) 
becomes 


u"(r)  +  Kb?r  =  0  ,  Vr  6  U  Qm  ■ 


The  displacement  field  can  be  determined  in  f 1m  and  Qm  up  to  a  constant 
(say,  cm  and  Cm )  by  solving  the  above  equation  and  using  the  boundary 
conditions.  This,  in  turn,  provides  an  expression  for  X\f(r)  and  A m(r)  from 
(23)1,3.  Furthermore,  solving  (23)2  gives  the  displacement  field  in  up  to  a 
constant  (say,  Ci).  With  this  displacement  field  one  can  solve  (24)  in  ft;  and 
determine  a(r)  up  to  a  constant  (say,  C2).  The  solution  is  displayed  more 
conveniently  in  nondimensional  form.  Define 


u  T 


a  =  —  ,  where  7  =  kui7!7  , 

a\[ 


and  the  following  nondimensional  parameters: 


_  2cr0 
;  “  aMKcjV7  ’ 


(25) 


The  parameter  /,  referred  to  as  the  loading  parameter ,  includes  the  applied 
loads,  the  upper  bound  of  the  material  properties  and  a  term  related  to  the 
body  force,  so  it  can  be  interpreted  as  the  ratio  of  surface  to  body  forces. 
The  parameter  /?,  referred  to  as  the  moduli  parameter ,  corresponds  to  the 
nondimensional  ratio  of  the  lower  and  upper  bounds  of  the  reinforcement  or, 
equivalently,  to  the  nondimensional  modulus  of  the  term  am  I  of  the  optimal 
material  Co-  The  displacement  field  and  the  maximum  principal  stiffness 
(reinforcement)  are  given  by 


r  |f(CM-f2)  2 

<  ^  (3C!  +  2f5) 

{  S  [cn  +  f  (3cq  -  r2)] 


f  6  ftM  =  [0,  f  1]  , 

f  6  ft;  =  (fi,f2J  , 
f  €  ftm  =  [f2, 1]  , 


u(f)  = 


(26) 
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and 

a(f) 


1 

C2 

\f?  ^3ci  +  2r  §  j 

P 


€  Qm  —  [0,  fi]  , 
€  O.i  =  [fi,f2]  i 

€  flm  =  (f2)  l]  , 


(27) 


where  the  constant  Co  in  (26)3  is  Co  =  y/l  +  f/P-  Two  relevant  nondimen- 
sional  quantities  can  be  identified  in  this  solution:  the  loading  and  moduli 
parameters  /  and  0  defined  in  (25)  which  can  be  used  to  characterize  the  dif¬ 
ferent  loading  cases.  The  constants  c™,  cm,  ci,  C2,  n  and  f2  are  determined  as 
follows:  since  fi  and  f2  represent  the  location  of  the  boundaries  between  Qm, 
Cli  and  ft,-,  respectively,  then,  for  an  admissible  solution,  it  is  required 
that  Aw(fi)  =  Am(f2)  =  0.  These  conditions  can  be  satisfied  by  choosing  cm 
and  Cm  as  follows: 

cM  =  (9  +  2 Vl5)f J  ,  ^  =  4^  (3co  -  l&Jff  +  7f\)  ■ 

Furthermore,  two  additional  relations  can  be  obtained  from  the  continuity  of 
u  at  f  =  f i  and  f  =  f2,  i.e., 


There  are  other  possible  values  for  cm  and  Ci  (cm  is  unique),  but  these  can  be 
discarded  a  posteriori  based  on  the  requirements  that  the  solution  should  not 
be  singular  and  that  the  multipliers  Am  and  A m  should  be  positive  in  [0,  rx) 
and  (f2,l]  respectively.  Finally,  the  value  of  c2  and  an  additional  relation 
between  fx  and  f2  can  be  obtained  by  enforcing  the  continuity  of  a  at  f  = 
and  f  =  f 2  which  gives 

c2  =  6(4  +  y/lb )t\12  =  fiy/fi  (1  +  -  2f2/2j  . 

The  values  of  f  i  and  f2  can  be  computed  numerically  from  the  above  relations. 
As  /  increases,  so  do  f\  and  f2.  For  some  values  of  /  ,  f2  >  1,  then  the  solution 
needs  to  be  modified  since  fim  =  0  (i.e.,  the  solution  never  reaches  its  lower 
bound  0).  In  that  case  the  solution  depends  on  the  loading  parameter  / 
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but  not  on  (5.  The  value  of  can  be  obtained  from  the  boundary  condition 
=  <7o,  which  can  be  expressed  as 

4(4  +  \/l5)ff  =/(2  +  (l  +  \/l5)f?)  .  (28) 

Finally,  for  /  large  enough,  fj  >  1,  hence  the  solution  is  trivial:  a  =  1  and 
the  boundary  condition  is  satisfied  with  cm  —  2>(f  +  1).  From  (27)  one  can 
obtain  the  total  resource,  Rt  =  a  J"j  adr,  where  a  is  the  cross-sectional  area. 
Since  a(f)  is  a  monotonically  decreasing  function  from  1  to  0  as  f  ranges 
from  0  to  1,  then  0  <  0  <  Rt  <  1  where  -Rt  =  Rj(oiMal).  Therefore, 
one  can  always  choose  an  upper  bound  R  =  R/(ctMal)  on  the  resource  such 
that  Rt  <  R  <  l,  i.e.,  such  that  the  resource  constraint  is  not  active  but 
Rt  ^  0.  This  confirms  the  assumption  that  there  are  non-trivial  solutions 
with  A  =  0.  Clearly,  if  the  prescribed  R  is  such  that  Rt  >  R,  then  (26)- 
(27)  would  not  be  an  admissible  solution  since  the  resource  constraint  would 
be  violated.  However,  an  analytical  expression  for  the  solution  when  the 
resource  constraint  is.  active  is  not  available. 

Figure  1  represents  the  total  resource  Rt  for  different  values  of  the  moduli 
parameter  0  and  loading  parameter  /.  Observe  that  for  values  of  /  large 
enough,  the  solution  does  not  depend  on  0  anymore  and  all  curves  merge 
into  a  common  envelope  which  can  be  thought  of  as  the  limit  case  0  -4  0  for 
all  values  of  /  (see  insert  in  Figure  1).  If,  for  a  given  pair  0,  /,  the  prescribed 
upper  bound  on  the  resource  constraint  R  is  above  the  corresponding  curve 
for  Ru  then  the  resource  constraint  is  not  active.  The  maximum  principal 
stiffness  a(f)  is  shown  in  Figure  2  for  different  values  of  0  and  a  common 
loading  parameter  /  =  0.1.  The  solid  curve  represents  the  case  for  which 
the  solution  does  not  depend  on  0  for  the  given  /  (i.e.,  from  Figure  1,  0 
is  below  some  critical  value  0.).  In  the  first  part  of  each  curve,  a  is  at  its 
upper  bound  then  it  decreases  monotonically  to  its  lower  bound  0  (except 
for  0  <  0.).  The  corresponding  minimum  values  of  the  objective  functional 
are  W.{0  =  0.3)  =  5.11  -lfT2,  W.{0  =  0.2)  =  4.07 -10"2  and  W.  =  3.30 -10"2 
for  0  <  0..  Although  0  is  viewed  as  a  fixed  parameter,  the  optimal  solution 
is  reached  for  values  of  0  <  0..  However,  as  mentioned  before,  0  needs  to 
be  bounded  away  from  zero  in  order  to  have  a  stable  structure.  Figure  3 
corresponds  to  the  displacement  field  as  a  function  of  f  for  different  values 
of  0.  Observe  that  at  f  =  1  the  displacement  of  the  solution  that  does 
not  depend  on  0  (solid  line)  is  greater  than  for  the  other  curves,  however, 
the  energy  norm  which  measures  structural  stiffness  is  smaller.  Similarly, 
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as  shown  in  Figure  4,  the  stress  distribution  is  also  smaller  for  ,3  <  /?., 
which  corresponds  to  the  stress-based  notion  of  optimal  structure.  Figure  5 
shows,  for  different  values  of  p,  the  (negative)  unconstrained  gradient  of  the 
Lagrangian  L  as  defined  in  (16)  (i.e.,  ~(L'  4-  Am  -  XM)h2  =  if  —  2fii, 
where  the  dot  represents  differentiation  with  respect  to  f).  For  each  curve, 
the  value  of  the  function  in  the  first  interval  Qm  =  [0,  f i)  (where  the  local 
upper  bound  is  active  (a  =  1)  and  where  L'  =  A,vf  =  0)  corresponds  to 
>  0;  in  the  second  interval  ft,-  =  [f1(f2]  (where  no  local  bounds  are 
active)  to  V  —  \m  =  Am  =  0  and  in  the  third  interval  Qm  =  (fr2, 1]  (where  the 
local  lower  bound  is  active  (<5  =  /?)  and  where  V  =  Am  =  0)  to  —  Am/7'  0- 

Observe  that  for  P  <  P.  (solid  curve),  there  is  no  third  interval  since  =  0. 

To  illustrate  that  this  procedure  can  be  used  in  a  more  general  setting,  a 
three-dimensional  example  was  solved  using  a  finite  element  program.  Con¬ 
sider  again  the  case  of  the  optimal  distribution  of  material  properties  in  order 
to  obtain  maximum  structural  stiffness  of  a  rotating  prismatic  structure  of 
length  l  and  square  cross-sectional  area  h2.  The  side  attached  to  the  axis  of 
rotation  is  modeled  as  a  clamped  end.  The  remaining  sides  of  the  structure 
are  subjected  to  the  following  loads: 

Case  1:  shearing  load: 

{  —  $  r°ez  f°r  T  =  ^  (Uniform  shear  stress) 

—  \  0  otherwise 

Case  2:  torsional  load: 

{-T0ez  for  r  =  l  ,  x  €  [-/i,  -h  4-  e)  , 

r0ez  for  r  =  /  ,  x  6  [h  —  e,  h]  , 

0  otherwise 

where  e  is  a  small  number,  the  origin  of  coordinates  is  in  the  mid-section 
and  x  is  measured  in  the  eg  direction.  The  numerical  values  used  here,  for 

illustration  purposes,  are  as  follows:  /  =  1  m,  h  —  0.1  m,  ecu  =  1010  Pa, 

am  =  10s  Pa,  k  —  10-7  s2  •  m-2  and  u  =  500  RPM.  In  Case  1,  To  =  106 
Pa  and  in  Case  2,  r0  =  1.5  •  106  Pa,  e  =  0.01  m.  In  both  cases  the  resource 
constraint  is  active  with  R  =  0.5lh2a\f  =  0.5  •  108  Pa-m3.  The  maximum 
principal  stiffness  a  =  a/aw  is  shown  for  the  loading  Case  1  (shearing  load) 
in  Figures  6  (without  a  body  force)  and  7  (with  a  body  force).  Figure  8 
represents  contour  plots  of  a  at  the  cross-section  f  =  0.75  for  the  load  Case 
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1;  the  section  on  the  left  corresponds  to  the  case  when  the  body  force  is 
absent  and  the  right  one  when  it  is  included.  Observe  that,  in  the  latter 
cross-section,  the  top  of  the  structure  is  in  compression,  thus  b  •  u  <  0, 
which  results  in  a  greater  reinforcement  than  the  bottom  which  is  in  tension 
(see  end  of  Section  3).  In  contrast,  when  the  body  force  is  absent,  the 
solution  is  symmetric  with  respect  to  the  r,  8  mid-plane,  as  shown  in  the 
left  cross-section  in  Figure  8  (compare  also  Figures  6  and  7).  The  heuristic 
interpretation  is  as  follows:  first,  observe  that  in  this  example  the  body  force 
always  creates  tension  at  any  point  of  the  structure.  Now,  consider  a  point 
of  the  structure  which  is  in  compression  because  of  the  applied  shear  load.  If 
the  material  is  locally  stiffer  (hence,  the  local  mass  density  is  higher),  then 
the  body  force  that  tends  to  oppose  compression  would  be  higher  resulting  in 
an  overall  smaller  displacement.  The  opposite  effect  occurs  at  points  which, 
because  of  the  applied  shear  load,  are  in  tension. 

For  the  loading  Case  2  (torsional  load),  Figure  9  shows  the  maximum 
principal  stiffness  when  the  body  force  is  not  present  and  Figure  10  when 
it  is.  In  Figure  11,  tyvo  cross-sections  at  f  =  0.75  are  shown;  the  left  one 
corresponds  to  the  case  without  the  body  force  and  the  right  one  when  the 
body  force  is  taken  into  account.  The  procedure  reproduces  the  well-known 
optimal  “circular"  layout  (within  the  limitations  of  the  mesh  and  the  effect  of 
the  boundary  conditions  at  r  =  0  and  r  =  /).  The  center  of  the  cross-section 
is  at  the  lower  bound  (hence,  it  can  interpreted  as  a  “weak”  region,  i.e.,  as 
a  “hole”).  This  illustrates  the  ability  of  the  method  to  obtain  the  optimal 
shape  as  a  by-product  of  the  material  optimization.  The  effect  of  the  body 
force  is  reflected  in  the  gradual  decrease  of  reinforcement  from  r  =  0  to  r  =  /, 
as  shown  by  the  cross-sections  in  Figure  11  (compare  also  Figures  9  and  10). 

5  Discussion  and  conclusions. 

The  method  presented  here  identifies  a  positive-definite  optimal  material  that 
maximizes  the  structural  stiffness  when  a  structure-dependent  body  force  is 
taken  into  account.  Two  interesting  aspects  arise  due  to  the  presence  of  the 
body  force:  as  opposed  to  problems  with  zero  body  force,  the  resource  con¬ 
straint  might  not  be  active  yet  it  is  possible  to  have  a  non-trivial  solution 
where  fn<p(a)dv  =  R.  ^  0.  This  suggests  that  there  is  a  “natural”  up¬ 
per  limit  R.  on  the  amount  of  reinforcement  material  that  needs  to  be  used 
(i.e.,  trying  to  specify  from  the  onset  a  greater  amount  will  not  affect  the 
optimal  solution).  Also,  it  is  not  necessarily  true  that  the  presence  of  the 
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body  force  will  result  in  a  local  increase  of  the  material  properties  in  order 
to  stiffen  the  structure.  In  fact,  as  shown  by  (21),  this  depends  on  whether 
the  scalar  product  of  the  body  force  and  the  displacement  is  negative  (local 
stiffening)  or  positive  (local  weakening).  It  is  usually  assumed  that  at  points 
on  the  surface  where  the  loads  are  applied,  the  optimal  layout  would  consist 
of  some  type  of  reinforcing  structure  to  support  those  loads.  However,  this 
is  not  the  case  if  the  structure  itself  turns  out  to  increase  substantially  the 
work  done  by  the  body  force.  Clearly,  as  shown  by  the  optimality  conditions, 
there  is  a  compromise  between  these  two  opposite  effects  and  the  outcome 
can  be  quantified  based  on  a  characteristic  magnitude  of  the  surface  loads, 
the  upper  bound  on  material  properties,  a  characteristic  magnitude  of  the 
body  force  and  whether  the  resource  constraint  is  active  or  not.  If  the  inertial 
forces  are  small  compared  to  the  applied  loads  (i.e.,  large  loading  parameter 
/)  then,  as  expected,  the  procedure  converges  to  a  solution  similar  to  the  case 
without  body  force.  However,  if  the  inertial  forces  are  large  (small  /),  then 
the  procedure  might  fail  to  develop  a  suitable  structure  to  support  the  loads 
and  large  displacement  gradients  could  occur  at  those  points.  Nevertheless, 
a  layout  which  includes  a  supporting  structure  for  the  loads  can  always  be 
obtained  by  selecting  an  appropriately  large  value  for  the  lower  bound  am. 
Finally,  it  is  worth  noting  that  the  numerical  implementation  of  this  method 
for  three-dimensional  problems  can  be  achieved  with  relative  ease  by  com¬ 
bining  (existing)  finite  element  programs  and  customized  optimizing  codes. 
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Rt=  Rtlctual 


f 

Figure  1:  One-dimensional  case:  normalized  total  resource  Rt  =  i?t/(ajv/Ql) 
for  different  values  of  the  loading  and  moduli  parameters  /  and  p.  The  insert 
corresponds  to  the  limit  case  p  — »  0  and  shows  the  full  range  0  <  Rt  <  1. 
For  prescribed  values  of  /  and  P ,  the  resource  constraint  is  not  active  if  the 
prescribed  R  is  above  the  corresponding  Rt  and  active  otherwise. 
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a  =  al  aM 


r 

Figure  2:  One-dimensional  case:  value  of  the  maximum  principal  stiffness 
a  (reinforcement)  in  the  radial  direction  for  various  values  of  the  moduli 
parameter  ft  and  a  common  loading  parameter  /  =  0.1.  The  solid  curve 
represents  the  limit  case  for  which  the  solution  does  not  depend  on 
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u  -  u/ly 


r 

Figure  3:  One-dimensional  case:  radial  displacement  u(f)  for  various  values 
of  the  moduli  parameter  /?  and  a  common  loading  parameter  /  =  0.1.  The 
solid  curve  represents  the  limit  case  for  which  the  solution  does  not  depend 
on  (3. 
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0  =  0/  aMy 
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Figure  4:  One-dimensional  case:  stress  distribution  as  a  function  of  radial 
position  for  various  values  of  the  moduli  parameter  (3  and  a  common  loading 
parameter  /  =  0.1  (i.e.,  the  boundary  condition  is  a(l)  =  // 2). 
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Figure  5:  One-dimensional  case:  (negative)  unconstrained  gradient  of  the 
Lagrangian  (=  u  —  2fu)  for  various  values  of  the  moduli  parameter  /?  and 
a  common  loading  parameter  /  =  0.1.  The  solid  curve  represents  the  limit 
case  for  which  the  solution  does  not  depend  on  /?.  See  text  for  detailed 
explanation. 
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Figure  6:  Loading  Case  1  (shearing  load):  contour  plot  of  the  maximum  prin¬ 
cipal  stiffness  (reinforcement)  without  a  body  force.  Observe  the  symmetry 
of  the  solution  with  respect  to  the  er,e$  mid-plane. 
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Figure  7:  Loading  Case  1  (shearing  load):  contour  plot  of  the  maximum 
principal  stillness  (reinforcement)  with  a  body  force.  The  upper  section  has 
a  greater  reinforcement  than  the  lower  section. 
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Shearing  load  r=  0.75 


No  rotation  Constant  rotation 


Figure  8:  Loading  Case  1  (shearing  load):  cross-sectional  contour  plot  of 
the  maximum  principal  stiffness  (reinforcement)  at  f  =  0.75;  the  left  plot 
corresponds  to  the  case  without  a  body  force  (symmetric  solution)  and  the 
right  one  to  the  case  with  body  force  (greater  reinforcement  in  upper  section 
which  is  in  compression,  lesser  reinforcement  in  lower  section  which  is  in 
tension). 
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Figure  9:  Loading  Case  2  (torsional  load):  contour  plot  of  the  maximum 
principal  stiffness  (reinforced  material)  without  a  body  force.  Observe  that 
the  distribution  of  reinforcement  is  prismatic,  except  close  to  f  =  0  and  f  =  1 
due  to  the  boundary  conditions. 
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Figure  10:  Loading  Case  2  (torsional  load):  contour  plot  of  the  maximum 
principal  stiffness  (reinforced  material)  with  a  body  force.  Observe  the  grad¬ 
ual  decrease  in  reinforcement  from  f  =  0  to  f  =  1 
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Figure  11:  Loading  Case  2  (torsional  load):  cross-sectional  contour  plot  of 
the  maximum  principal  stiffness  (reinforcement)  at  f  =  0.75;  the  left  plot 
corresponds  to  the  case  without  a  body  force  and  the  right  one  to  the  case 
with  body  force. 
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Abstract 

An  effective  procedure  for  the  solution  of  the  optimal  layout  and  properties  of  an  elastic 
material  is  presented.  Design  of  2D  and  3D  structures  is  considered  for  minimum 
compliance  under  a  single  load  condition  with  a  constraint  on  the  amount  of  material  and 
bounds  on  material  properties.  The  method  is  based  on  an  already  established  formulation 
for  continuum  design  with  the  elasticity  tensor  being  the  design  variable.  The 
optimization  with  respect  to  local  dependence  of  the  tensor  on  orientation  may  be 
separated  from  the  variation  of  material  properties  density  (represented  by  an  invariant  of 
the  tensor)  over  the  field  and  solved  analytically.  A  weighting  function  is  introduced  to 
measure  the  cost  of  the  material  properties  density.  An  optimal  layout  composed  of  either 
zero  or  full  material  properties  density  (optimal  topology)  is  predicted  out  of  a  sequence 
of  solutions  to  materia]  properties  design  sub-problems,  each  corresponding  to  a  different 
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weighting  function.  As  a  first  step,  the  optimal  distribution  of  continuously  varying 
material  properties  density  is  predicted  for  uniform  weighting  function.  For  each 
subsequent  sub-problem,  the  weighting  function  is  systematically  adjusted  and  applied  to 
the  optimal  material  distribution  obtained  in  the  previous  step.  This  sequence  results  in  a 
segregation  of  the  material  properties  density  to  either  its  lower  or  upper  bound.  The  final 
design  can  be  identified  as  topology  design,  where  the  lower  bound  value  is  small 
enough — on  the  order  of  10'5  in  the  examples — compared  to  the  value  for  the  upper 
bound  (full  material).  Thus,  this  technique  consists  in  a  gradual  transition  from  the  initial 
optimal  continuous  elastic  structure  to  one  having  sharply  defined  topology  by  the  use  of 
a  weighted  total  resource  constraint.  The  method  was  successfully  implemented  using  a 
commercial  finite  element  analysis  package  to  solve  several  example  problems  in  2D  and 
3D. 


1.  Introduction 

This  study  centers  on  the  development  of  a  procedure  for  the  prediction  of  the  optimal 
material  distribution  and  properties  of  two-  and  three-dimensional  linear  elastic 
continuum  structures  using  finite  element  discretization.  The  goal  is  to  find  the  least 
compliant  design  subject  to  a  constraint  on  the  total  amount  of  available  material.  The 
material  layout  of  the  continuum  structure  design  is  described  by  a  density  field  that 
assigns  value  to  a  pointwise  measure  of  the  material  properties.  Restricting  the  material 
properties  density  function  to  take  one  of  two  discrete  values  makes  it  possible  to  treat 
the  topology  optimization  problem.  The  upper  bound  of  the  material  density  is  identified 
as  the  fully  concentrated  measure  of  the  structural  material,  while  the  lower  bound 
represents  extremely  weak  material — essentially  a  lack  of  structure.  Thus  the  contour  of 
the  concentrated  material  properties  density  describes  the  topology  of  the  optimized 
structure. 

Most  of  the  well-established  characterizations  for  the  topology  design  within 
structural  optimization  are  based  on  the  assumption  that  the  material  properties  are 
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represented  by  a  composite  composed  of  two  isotropic  elastic  components.  Material 
layout  and  properties  optimization  as  well  as  topology  optimization  problems  have  been 
treated  using  such  formulations.  It  follows  from  the  requirement  of  well-posedness  and 
numerical  tractability  that  models  of  this  type  have  to  be  relaxed,  since  the  optimal  design 
is  often  composed  of  an  infinitely  fine  composite  microstructure.  Choosing  the  procedure 
to  represent  effective  elastic  moduli  of  possible  material  microstructures  is  the  usual  way 
to  achieve  relaxation.  In  some  schemes  all  possible  arrangements  are  present  (see  e.g. 
Diaz  &  Lipton  (1997)),  while  in  others  only  a  subclass  of  microstructures  is  allowed  as  a 
consequence  of  specific  composite  microgeometry  parametrization.  Regardless  of 
possible  microstructure  restrictions,  relaxation  is  usually  achieved  using  a 
homogenization  technique.  More  details  on  the  procedures  and  models  for  topology 
design  related  to  the  two-material  composite  approach  can  be  found  in  Bendsde’s  treatise 
(1995)  and  in  the  proceedings  Bendspe  and  Mota  Soares  (1993). 

In  this  paper  we  use  the  formulation  of  structural  optimization  where  the  design  is 
represented  by  the  unrestricted  elastic  properties  tensor  field.  There  is  no  restriction  on 
the  elasticity  tensor  except  for  those  required  to  represent  stable  physical  materials  and  to 
provide  existence  of  elastostatics  problem  solutions.  This  formulation  was  originally 
presented  in  Bendspe  et  al.  (1994)  for  linear  materials  and  subsequently  extended  for 
multiobjective  optimization  problems  (Bendspe  et  al.  (1995)),  nonlinear  materials 
(Bendspe  et  al.  (1996))  and  optimal  design  of  a  structure  with  design-dependent  field 
force  (Turteltaub  and  Washabaugh  (1998)).  It  is  noteworthy  that  optimization  of  the 
material  tensor  orientation  forms  the  local  problem  and  may  be  separated  from  the 
determination  of  the  material  properties  density  distribution  over  the  structure.  The  local 
problem  has  an  analytical  solution  (Bendspe  et  al.  (1994)).  Thus  the  original  problem  is 
effectively  reduced  to  one  of  finding  the  optimal  field  of  an  invariant  of  the  elasticity 
tensor  (the  trace  of  the  elasticity  tensor  is  used  here  for  this  purpose).  The  resulting 
material  density  distribution  is  continuous  in  the  design  domain  and  non-uniform  except 
for  those  regions  where  local  bounds  are  reached.  Note  that  subregions  bounded  by 
adjacent  contours  on  a  plot  of  results  for  the  density,  where  each  contour  corresponds  to  a 
level  of  density,  represent  a  form  of  design  topology.  This  property  is  exploited  for 
prediction  of  optimal  structural  topology  in  the  form  of  a  sharp  (black/white)  image 
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(Guedes  and  Taylor  (1997)).  There,  the  usual  form  for  the  total  resource  constraint  is 
replaced  by  a  generalized  form,  having  a  weighting  function  introduced  to  measure  the 
local  cost  of  material  properties  density  (trace)  of  the  elasticity  tensor.  (A  generalized 
formulation  of  unit  cost  fields  is  described  in  Taylor  and  Washabaugh  (1995)  and 
(1997).)  Applying  stepwise  adjustment  of  the  weighting  leads  to  a  reduction  in  regions 
having  intermediate  density  levels  in  redesign  results  for  optimal  material  density.  Only 
two  values  of  material  properties  density  remain  in  the  refined  topology  final  design. 
One,  identified  by  black,  corresponds  to  full  measure  of  the  material  while  the  other 
(white)  has  the  trace  measure  at  the  lower  bound  and  symbolizes  lack  of  material. 

The  free  material  tensor  design  problem  is  well  posed  and  does  not  require 
relaxation  as  is  needed  in  formulations  based  on  two-material  composite  models.  The 
simplicity  of  the  material  tensor  approach  in  the  optimal  material  properties  field  problem 
is  preserved  for  the  generalized  formulation  with  weighted  resource  constraint  as  used  for 
sharp  image  (black/white)  design.  Consequently  the  usual  steps  required  in  existing 
procedures  for  topology  design,  such  as  imposition  of  additional  constraints  on  the  design 
space,  or  penalization  techniques  for  design  variable,  are  avoided  here. 

In  this  paper  we  deal  with  the  sharp  image  design  using  the  characterization 
which  closely  follows  that  described  by  Guedes  and  Taylor  (1997).  A  computational 
procedure  based  on  the  elastic  free  material  tensor  approach  for  prediction  of  optimal 
material  layout  (topology)  was  implemented  using  a  commercial  finite  element  analysis 
package  and  tested  on  2D  and  3D  examples.  Topology  design  results  for  several 
problems  are  presented  in  this  paper.  Several  objectives  represent  the  core  of  the  present 
work:  it  is  shown  that  the  method  for  predicting  optimal  topology  can  be  extended  to  a 
three-dimensional  setting  with  virtually  no  modifications.  Furthermore,  the  examples 
included  in  the  present  paper  illustrate  situations  where  the  weighted  resource  constraint 
method  can  be  effectively  applied  in  order  to  satisfy  geometrical  requirements  of  a  design 
(e.g.,  formation  of  holes  in  predetermined  regions).  Of  practical  relevance,  it  is  also 
shown  that  the  method  can  be  easily  implemented  using  available  finite  element 
programs. 
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2.  Problem  statement 


The  analytical  model  for  the  optimal  design  of  linear  structures  with  an  isoperimetric 
(cost)  constraint  presented  here  closely  follows  the  formulation  given  by  Bendspe,  et  al. 
(1994)  and  subsequently  elaborated  to  predict  optimal  topology  by  Guedes  and  Taylor 
(1997).  It  differs  however,  from  the  original  by  means  of  how  the  specified  unit  cost  field 
is  defined.  Here  the  unit  cost  is  required  to  belong  to  the  set  of  functions  such  that  the 
product  of  a  function  and  material  properties  density  is  normalized  to  the  total  material 
resource.  With  this  restriction,  the  original  total  resource  constraint  is  automatically 
satisfied  whenever  the  constraint  of  total  design  cost  is  satisfied. 

We  consider  a  linear  elastic  continuum  structure  occupying  a  domain  Cl,  subjected 
to  the  single  set  of  body  forces  /  ,  boundary  traction  t  (on  part  of  the  boundary  dClt ) 
and  zero  displacement  on  c) Clu .  The  design  variables  are  the  local  elastic  properties, 
which  are  described  by  the  material  modulus  tensor  E .  An  invariant  measure  (the  trace 
p  of  the  modulus  tensor,  i.e.,  p  =  E~iJ)  is  subject  to  local  and  global  constraints  (see 
(lc,d)  below).  We  seek  to  predict  the  structure  (in.  terms  of  its  material  properties 
distribution)  that  satisfies  equilibrium  and  minimizes  the  total  work  of  external  forces  at 
this  state  (i.e.,  minimum  compliance  design).  This  problem  can  be  expressed  as 


[PI] 


nun 


1(“) s  J  fiuj&Cl  + 1  ► 

n  <n, 


subject  to: 

|^(u)£(/u£u(v)da  =  /(v)  Vve  t/, 

n 


E>  0, 

|  wpdCl  -  R  <  0 , 
n 

p<p<p. 


(1) 

(la) 

(lb) 

(lc) 

(ld) 


Here  u,v  represent  displacement  fields  belonging  to  a  set  U  of  kinematically  admissible 
fields,  E  or  symbolizes  the  fourth-order  symmetric  positive  semi-definite  elasticity 
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tensor  field  and  £-(«)  =  ^(m; .  +«;j  )  stands  for  the  strain  tensor  (associated  with  either  u 
or  v).  In  the  third  and  fourth  constraints  (lc,d),  w=w(x)  >  0  is  a  given 
function — referred  to  as  the  weight  function — whereas  p  and  p  are  specified  local 
bounds  on  the  trace  p .  The  global  cost  of  the  design  is  limited  by  a  specified  value  R . 

Optimal  material  properties. 


The  equilibrium  state  described  by  equation  (la)  in  weak  form  can  be  expressed 
alternatively  as  the  minimum  of  the  potential  energy  with  respect  to  admissible 
displacement  fields  u.  Furthermore,  minimization  of  the  total  work  of  external  forces 
stated  in  (1)  might  be  exchanged  by  maximization  of  the  potential  energy  because  the 
total  work  of  external  forces  at  equilibrium  equals  the  negative  of  twice  the  potential 
energy.  Thus  the  original  problem  [PI]  can  be  reformulated  as 

(2) 


max  min 

£  «i/ 


subject  to 


E>  0, 

|  wpdCl  —  R<  0 , 
n 

P^p<p, 


(2a) 

(2b) 

(2c) 


where  the  potential  energy  is  given  by 


n(£,H)={j£(f(ii)£#ufH(H)dn  -/(«), 


(3) 


(See  also  Bendspe  et  al.  (1994).)  The  maximization  with  respect  to  the  design  variable 
(i.e.,  the  elasticity  tensor  E )  can  be  separated  into  two  sets  of  parameters:  one  that 
corresponds  to  the  local  material  properties  orientation  and  another  which  refers  to  the 
variation  over  the  field.  The  latter,  representing  the  cost  of  the  material,  is  taken  to  be  the 
first  invariant  p  of  the  elasticity  tensor  as  in  the  global  cost  constraint  (2b). 
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After  this  separation  of  the  design  maximization  into  two  parts,  our  problem  takes 
the  form 


max 

p 

subject  to :  p<  p<p 

J  wpdQ  -  R  <  0 

n 

Interchanging  the  inner  max  and  min  problems  (which  is  valid  since  it  is  a  saddle 
point  as  proved  in  Bendspe  et  al.  (1994))  makes  it  possible  to  recognize  the  optimal  local 
material  properties  in  the  form 


max 

£ 


min  n  (p,E,u). 

1*6  U 


(4) 


subject  to:  E  >  0 

£....  =  p 
w  r 


Em=P 


(5) 


This  corresponds  to  an  orthotropic  material  having  axes  of  orthotropy  co-aligned  with  the 
principal  strain  directions.  According  to  (5),  the  strain  energy  stored  in  the  optimized 
material  is  given  by 


U=lEijtJ^k!  =  T  P£ij£ij  •  (6) 

Thus  the  optimal  material  has  the  same  effective  energy  of  an  isotropic  linear  elastic 
material  with  zero  Poisson  ratio  and  Young  modulus  equal  to  p .  This  feature  is  a 
consequence  of  the  specific  orthotropy  orientation  (i.e.,  it  is  related  to  the  coupling 
between  the  elastic  properties  and  the  strain  field  at  equilibrium).  However,  the  cost 
measured  by  the  trace  of  the  elastic  moduli  tensor  of  the  optimal  material  is  six  times 
lower  than  for  the  isotropic  one  (i.e.,  the  trace  of  p{6ik6jl  +  is  6p ) 

The  analytical  result  for  the  local  optimal  material  properties  (5),  (6)  can  be 
incorporated  into  the  problem  statement  (4).  The  number  of  design  variables  is  then 
reduced  to  only  one,  which  is  the  distribution  of  the  trace  of  the  elasticity  tensor  p  .  Thus 
the  problem  statement  may  be  expressed  as 
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[P2] 


max  min-! 

p  ueU 


t | - J /;•  utdn -  \l «,ds 

n  n  ao, 


subject  to: 
p-p<0, 
p~p<0. 


(7) 

(7a) 

(7b) 


Jwpdn-/?<0.  (7c) 

a 


According  to  (7),  the  original  problem  for  the  prediction  of  the  optimal  material 
properties  field  consists  of  obtaining  the  optimal  distribution  of  the  trace  of  the  elasticity 
tensor  and  the  corresponding  strain  field.  Applying  the  analytical  solution  (5)  of  the  local 
problem  gives  the  pointwise  optimal  orientation  of  the  elasticity  tensor.  We  note  again 
that  the  original  problem  (see  Bendsde  et  al.  (1994))  can  be  obtained  from  the  present 
formulation  simply  by  setting  w(r)=  1 . 

Optimality  conditions. 


The  Lagrangian  corresponding  to  the  constrained  design  problem  (P2)  is  given  by 


L  =  I J  pe&jdO.  -/(«)-  J  A  (p  -  p  )d£2  -  J  X{p  -  p) dO 

n  on 

JwpdO-i?"! 


(8) 


and  the  optimality  condition  for  the  distribution  of  elastic  material  properties  p  may  be 
stated  as 


^  -  X  +  X  -  Aw  =  0 . 

The  solution  also  satisfies 


(9) 
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(10) 


A  >0  A(p- p)  -  0, 

A>0  A(  p-p)  =  0. 


The  resource  constraint  is  satisfied  as  an  equality  (see  Bends0e  et  al.  (1994))  and  the 
corresponding  multiplier  is  positive,  i.e., 

A>0  and  Jwpdn-tf  =  0.  (11) 

n 


Weighting  function. 

Henceforth,  we  consider  a  specific  stepwise  constant  weighting  function,  defined  as 
follows: 

w(x)  =  (vv  -  w)  l(p(x),  pj  +  w  (12) 


where  l{p(x),pth)  is  the  unit  Heaviside  function,  i.e., 


fo 


l(p(x),pth)  =  \ 


1 


for  p(x)  <  plh 
for  p(x)  >  plh 


(13) 


and  plh  is  a  designated  cutoff  value  between  p  and  p  .  From  (9)  and  (10),  the  optimality 
condition  in  the  subdomain  where  p  <  p(x)<J)  becomes 


\£..£;j=hw  Vxg  Q  s.t.  p(x)e  (£,p)  , 


(14) 
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hence  ^  £()£(y  =  A  for  p  <  p  <  p,A  and  =  A  for  pth<  p<  p  .  This  condition 

states  that  the  weighted  unit  strain  energy  for  the  optimal  material  configuration  should 
be  uniform  in  the  design  domain  Q  except  possibly  for  intervals  or  points  in  Cl  where  p 
equals  p  or  p  .  The  strain  energy  is  to  be  taken  at  the  equilibrium  state  (inner  max  in 

P2).  A  jump  of  the  unit  strain  energy  i  £„£„  occurs  at  the  boundary  between  low  and  high 

weighted  regions  as  the  consequence  of  the  discontinuous  specific  weighting  function 

(12). 

3.  Computational  procedure 

The  optimal  material  properties  design  problem  (P2)  was  implemented  and  solved  for  2D 
and  3D  examples  using  a  finite  element  approximation  for  a  discrete  representation  of  the 
continuum.  The  optimal  material  properties  are  given  by  the  strain  field  at  equilibrium 
(5),  but  in  practice  it  is  simpler  to  use  an  energetically  equivalent  isotropic  material  (see 
equivalent  form  of  strain  energies  (6)).  Standard  procedures  supported  by  a  FEM 
analysis  package  were  used  to  obtain  the  displacement  field  at  equilibrium  (internal  min 
in  (P2)).  To  solve  the  external  max  problem,  an  optimality  criterion  method,  based  on 
conditions  (6)-(9),  was  applied.  In  order  to  predict  the  optimal  topology  of  the  structure,  a 
stepwise  procedure  was  developed.  This  procedure  consists  on  solving  a  sequence  of 
optimization  problems  (P2)  for  a  given  weighting  function,  which  is  updated  at  the  end  of 
each  solution.  The  procedure  is  summarized  as  follows. 

As  a  first  step — to  initialize  the  topology  prediction — the  optimal  distribution  of 
elastic  material  properties  p0  is  computed  for  a  uniform  weighting  function  (i.e.,  with 
w0(r)=l).  Then,  a  cutoff  value  of  material  properties  close  to  p  is  selected  (i.e., 
p,°  =  p  +  A p,° ,  where  A p,°  is  a  small  positive  increment).  A  new  weighting  function 
w,(x)  =  w,(p0(:c))  is  obtained  from  (12)  and  (13)  based  on  the  optimal  solution  p0  and 
the  cutoff  value  p°  .  Next,  the  optimization  problem  for  this  new  weighting  function  w, 
is  solved,  resulting  in  a  new  optimal  distribution  of  material  properties  p, .  The 
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distribution  px  is  characterized  by  a  “shift”  in  p :  the  values  of  px  are  lower  in  regions 
where  initially  p0  <  p°h  and  higher  where  p0  >  pt°  .  In  subsequent  steps  the  cutoff  value 

is  increased,  a  new  weight  function  is  assigned,  and  the  optimal  design  problem  is  solved 
again.  This  scheme  is  repeated  until  only  the  lower  and  upper  bound  values  on  material 
properties  density  remain  in  the  final  design.  (The  computational  procedure  terminates 
when  the  number  of  elements  with  p  laying  in-between  bounds  is  lower  than  a  user- 
specified  number.)  In  other  words,  the  solution  procedure  results  in  a  gradual 
transformation  of  the  initial  optimal  continuous  structure  to  a  discrete  one,  the  latter 
obtained  via  a  sequence  of  optimization  subproblems.  The  solution  of  each  subproblem 
becomes  the  initial  design  for  the  next  optimization  problem.  The  entire  computational 
procedure  is  described  in  the  form  of  a  flow  diagram  (Fig.  1).  In  the  diagram,  the  index  n 
refers  to  the  iterative  procedure  to  find  an  optimal  distribution  of  material  properties  for  a 
given  weight  function  (i.e.,  problem  (P2)).  The  index  i  refers  to  an  internal  loop  to  satisfy 
the  global  resource  constraint  (7c).  The  optimality  criterion  method  used  to  solve 
problem  (P2)  involves  a  local  step  size  parameter  for  design  variable  change,  denoted  by 
a  (a  value  between  0  and  1).  The  convergence  tolerance  for  the  objective  function 
(compliance  /(«))  is  a  given  value  8,  whereas  the  tolerance  for  the  resource  constraint  is 
p.  Finally,  the  index  /  refers  to  the  sequence  of  optimization  problems  (P2).  After  each 
optimization  problem  l  (for  a  fixed  weighting  function),  the  new  cutoff  value  is  obtained 
as  p]h  =  pip  +  A p'lh ,  where  A p\h  is  the  cutoff  value  increment. 

Two  parameters  strongly  affect  the  procedure  and  its  result:  the  relative  weight 
ratio  wR  =  w  /  w  and  the  cutoff  value  increment  Apth .  The  first  of  them  determines  the 
dispersion  of  material  properties  p  in  a  single  step,  while  the  other  controls  how  fast  the 
weighting  function  changes  from  step  /  to  step  /  + 1 .  Numerical  experiments  showed  that 
usually  for  a  given  problem  there  exists  a  minimum  value  of  the  relative  weight  ratio  vvfi 
that  is  satisfactory  to  obtain  a  discrete  black  and  white  design  for  an  appropriate  cutoff 
value  increment.  For  weight  ratios  lower  than  wR ,  a  significant  number  of  elements 

where  p  lies  between  the  bounds  may  be  present  in  the  final  design,  regardless  of  A plh 
(i.e.,  wR  needs  to  be  sufficiently  large  in  order  to  reach  a  black/white  design).  On  the 
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other  hand,  there  is  no  unique  criterion  for  the  selection  of  the  cutoff  value  increment 
A pth .  The  sensitivity  of  the  design  with  respect  to  Apth  depends  strongly  on  the  cutoff 
value  pih  at  which  this  increment  is  applied.  Therefore,  the  rate  of  adjustment  of  the 
weighting  function  is  one  of  the  keypoints  of  the  algorithm  since  it  affects  the  quality  of 
the  final  design  and  the  computational  time.  The  method  adopted  here  is  as  follows:  let 
Vs(p'h)  be  the  volume  of  subdomain  for  which  the  trace  pt  at  step  /  is  below  the  /-th 
cutoff  value,  i.e., 

Vs(p'J=  jw  0,  =  (*)</>:}  .  (15) 

Recall  that,  from  (12)  and  (13),  Qfi,  essentially  correponds  to  the  region  where  the 

weighting  function  is  equal  to  its  lower  value  w.  To  adjust  the  cutoff  value  for  the  next 
optimization  subproblem  (step  l  + 1 ),  the  basic  idea  is  assume  a  value  for  Ap'lh  and  check 
that  the  volume  Vs  (p'lh  +  Ap1^ )  is  bounded  above  and  below  as  follows: 

J  (16) 

"  Vs(p)-Vs(p'lh) 

where  £<£,  £,£  e  (0,l),  are  user-defined  bounds.  The  cutoff  value  increment  is 
calculated  iteratively  in  order  to  satisfy  (16).  In  (16),  Vs(p )  refers  to  the  volume  of  the 
domain  occupied  by  the  material  that  is  not  at  the  upper  bound  p  .  The  advantage  of  this 
volume  criterion  (to  adjust  the  cutoff  value)  is  that  it  is  based  on  an  a  priori  known 
quantity  directly  related  to  the  topology  of  the  final  design. 


4.  Examples 

A  finite  element  discretization  was  used  for  the  numerical  solution  of  2D  and  3D 
problems  of  optima]  material  properties  layout  with  a  given  weighting  function.  The 
finite  element  analysis  package  ANSYS  was  incorporated  in  the  computational  procedure 
described  in  the  previous  section.  Standard  linear,  4-node  quadratic  elements  for  plane 
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problems  and  8-node  brick  elements  for  a  three-dimensional  example  were  applied.  The 
material  properties  may  vary  from  element  to  element,  but  they  are  assumed  to  be 
constant  inside  each  element. 

Example  1 

The  objective  in  the  first  example  is  to  design  the  stiffest  plane  structure.  The  structure  is 

restricted  to  a  rectangular  domain  of  medium  aspect  ratio  (3x5)  (see  Fig.  2).  A  state  of 

plane  stress  is  assumed.  Boundary  conditions  are  specified  at  the  left  edge  by  zero 

vertical  and  horizontal  displacements.  On  the  opposite  side,  a  concentrated  load  acts  in 

the  middle  of  the  edge.  The  displacement  field  is  approximated  by  a  quadrilateral  4-node 

element  and  a  100x60  uniform  mesh.  The  material  properties  are  constant  within  each 

element.  The  upper  bound  of  the  trace  of  the  elasticity  tensor  is  assumed  to  be  105  times 

higher  than  the  lower  bound  (i.e.,  wR  =  105 ).  The  volume  fraction  of  the  full  material  at 

the  final  design  is  specified  to  be  40%.  Thus  the  total  resource  of  material  properties  is 

R=QAp\dV.  ■ 
a 

Figure  la  shows  the  contour  plot  of  the  trace  of  the  elasticity  tensor  p  for  the  initial 
optimization  problem  with  a  uniform  weighting  function  (i.e.,w(r)=  1 ).  The  black 
contours  correspond  to  full  material  properties  (i.e.,  p-'p).  The  lighter  shades  of  gray 
indicate  lower  values  of  the  elastic  properties  trace.  Compared  to  the  initial  uniform 
material  distribution,  the  total  compliance  was  decreased  by  1/3  in  31  iterations  after 
which  convergence  of  10^  was  reached  (10'3  after  5  iterations).  Figures  lb  and  lc  show 

the  evolution  of  the  continuous  design  towards  black  and  white  for  slow  ( £  =  0.02 )  and 

rapid  (£  =  0.05)  cutoff  value  adjustments.  In  the  first  case,  the  finer  structure  features  a 
22%  increase  in  the  objective  function  (compliance)  relative  to  the  optimal  continuously 
varying  material  distribution  design.  (It  is  worth  recalling  that  the  black/white  solution  is 
sub-optimal  compared  to  an  optimal  structure  with  continuously  varying  properties.)  The 
second  one  is  28%  less  stiff  than  the  continuous  design.  We  note  that  despite  using  linear 
4-nodes  elements,  no  checkerboard  pattern  has  occurred.  However,  for  extremely  low 
cutoff  value  adjustment,  this  feature  was  observed. 
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Example  2 

The  problem  of  a  circular  plate  with  a  central  hole  compressed  along  a  diagonal  is  now 
considered.  The  ratio  of  the  external  to  the  internal  radius  equals  to  2.  The  load  is 
distributed  according  to  a  sinusoidal  function  on  1/10  of  the  external  edge.  A  state  of 
plane  stress  is  assumed.  The  structure  and  the  loading  have  two  axes  of  symmetry;  hence 
the  design  domain  can  be  restricted  to  one  quarter  of  the  original  domain  (see  Fig.  3). 
Again,  the  4-node  elements  were  used  to  discretize  V*  of  the  structure.  The  mesh 
consisted  of  5000  elements.  As  before,  the  same  local  bounds  on  the  trace  of  the  elasticity 
tensor  were  applied.  To  show  the  influence  of  the  total  resource  available  on  the  final 
result,  four  values  corresponding  to  20%,  30%,  40%  and  50%  of  the  volume  fraction 
were  applied  (keeping  the  local  constraints  unchanged).  Results  are  shown  on  Fig.  3.  The 
values  of  the  objective  function  for  the  optimal  designs  with  a  uniform  cost  function  (i.e., 
w(r)=  1 )  and  increasing  amount  of  available  material  (4  volume  fractions  stated  above) 

are,  in  each  case,  0.55,  0.58,  0.62  and  0.66  of  the  initial  values  corresponding  to  a 
uniform  material  distribution  (i.e.,  p  =  constant).  The  compliance  for  black/white  designs 

decreases  respectively  to  0.67,  0.70,  0.71  and  0.70  of  the  initial  values  for  uniform 
material  distribution. 

The  interesting  feature  observed  here  is  that  the  light  truss  structure  obtained  for 
the  least  available  material  case  (20%  of  volume  fraction)  evolves  to  a  solid  design  for 
increasing  total  resource  of  material. 

Example  3 

This  example  is  based  on  Lame’s  problem  in  plane  stress.  The  same  design  domain,  as  in 
example  2,  is  given.  The  load  acting  on  the  external  edge  is  a  uniform  pressure.  Thus  the 
problem  is  axisymmetric  for  continuous  distribution  of  material  properties  (i.e.,  it  is 
essentially  a  ID  problem  since  it  has  radial  symmetry).  However,  due  to  the  specific 
choice  of  a  weight  function  in  this  example,  we  expected  a  non-axisymmetric  solution  for 
the  discrete  design.  Therefore,  a  quarter  of  the  design  domain  was  chosen  for  this 
analysis.  The  structure  was  discretized  using  quadrilateral  4-node  elements.  The  total 
resource  constraint  was  specified  to  satisfy  40%  of  the  volume  fraction.  The  solution  for 
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the  uniform  weight  function  (i.e.,  w(x)=  1)  features  a  relatively  uniform  material 
properties  distribution  (see  Fig.  4).  The  elements  on  the  inner  edge,  having  full  material 
properties,  form  a  strong  reinforcing  ring.  However,  the  lower  bound  of  the  elasticity 
tensor  trace  is  not  reached  anywhere  in  the  design  domain. 

For  the  topology  design  problem,  the  procedure  was  very  unstable  for  different 
choices  of  the  design  parameters.  In  this  case,  the  uniform  material  properties  distribution 
was  apparently  the  reason  that  it  failed  to  predict  a  discrete  topology.  One  of  the 
pathologies  experienced  was  an  axisymmetric  checkerboard  pattern.  To  avoid  difficulties 
caused  by  the  uniform  material  distribution,  an  initial  user-defined  weighting  function 
was  designated.  This  function,  as  shown  in  Fig.  4b,  has  three  sectors  where  w  =  w  (and 
w  =  w  elsewhere).  The  final  topology  design  was  obtained  using  the  same  procedure  as 
in  the  previous  examples.  Even  though  this  topology  depends  on  the  initial  choice  of  w, 
the  introduction  of  arbitrary  weighting  functions  in  the  design  procedure  may  be  used  in 
practice  to  set  preferences  of  the  final  topology.  As  shown  in  this  example,  if  holes  are 
required  in  each  of  the  three  sectors  shown  in  Fig.  4b,  this  can  be  achieved  simply  by 
specifying  a  suitable  weighting  function. 

Example  4 

This  example  consists  of  a  3D-cantilever  beam  with  a  transverse  force,  similar  to  the  2D 
beam  shown  in  example  1.  The  structure  is  required  to  support  a  concentrated  load  at  the 
middle  of  the  face  opposite  to  the  clamped  face  (see  Fig.  5).  The  design  domain  has  the 
shape  of  a  3x3x5  brick  and  the  material  is  required  to  occupy  at  most  40%  of  its  volume. 
The  optimal  topology  is  shown  in  Fig.  5  (the  final  design  is  shown  in  gray  to  see  the 
mesh,  however  all  visible  elements  are  at  the  upper  bound  of  the  trace  of  the  elasticity 
tensor). 

Unlike  its  2D  analog,  the  resulting  structure  consists  of  two  parallel  “walls” 
separated  by  a  long  hole  inside  the  beam.  The  objective  function  is  35.8%  lower  than  the 
one  given  by  a  structure  of  uniformly  distributed  material  properties.  By  comparison,  the 
optimal  structure  with  continuously  varying  material  properties  and  same  volume  fraction 
provides  a  38.5%  reduction. 
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The  same  problem  was  subsequently  solved  for  volume  fractions  of  20%  and 
10%.  The  resulting  structure  for  a  volume  fraction  of  20%  is  topologically  similar  to  the 
one  described  above  (see  Fig.  6)  and  gives  a  40%  reduction  in  the  objective  function 
(compared  to  a  47%  reduction  for  the  structure  of  continuously  distributed  properties). 
For  10%,  however,  the  supporting  walls  were  broken  and  a  structure  close  to  the  coarse 
2D  design  was  obtained  (Fig.  7).  The  objective  function  in  this  case  is  much  lower  for  the 
continuous  structure  (48.5%)  than  for  the  black/white  structure  (35%). 

To  test  the  capacity  of  a  “common”  small  computer  system  to  solve  a  relatively 
large  scale  problem,  the  3D  example  was  solved  using  20-node  brick  elements  (instead  of 
8-node  bricks)  which  resulted  in  50  000  degrees  of  freedom  for  the  FE  model.  The 
solution  took  approximately  20  hours  on  a  Pentium  150  PC  machine.  The  results  were 
topologically  similar  to  those  discussed  before  and  the  difference  in  objective  function 
was  less  than  2%. 

It  has  been  observed  that  the  computational  procedure  shows  a  small  tendency  to 
develop  an  unstable  checkerboard  pattern.  However,  checkerboarding  occurs  only  for 
extremely  low  cutoff  value  adjustments.  This  feature. was  specifically  visible  when  the 
discretization  involved  linear  finite  elements  with  constant  material  properties  and  can  be 
eliminated  by  a  more  accurate  approximation  of  the  displacement  field  and  the  material 
properties  distribution  (Jog  et  al.  (1993)). 

5.  Conclusions. 

A  method  based  on  the  arbitrary  (free)  elastic  tensor  formulation  for  the  minimum 
compliance  problem  is  used  here  to  predict  the  optimal  material  layout  and  properties  in 
the  sense  of  structural  topology  design.  The  transformation  of  the  continuously  varying 
optimal  material  properties  field  to  a  discrete  layout  design  is  achieved  by  solving  a  series 
of  optimization  problems  with  weighted  total  resource  (cost)  constraint.  Following  this 
approach,  the  computational  procedure  was  developed  and  implemented  in  the  finite 
element  analysis  package  ANSYS.  The  results  of  calculations  performed  on  several  two 
and  three-dimensional  examples  showed  that  the  present  method  is  capable  to  solve 
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medium  and  large  problems,  being  relatively  more  efficient  than  familiar  schemes  for 
topology  design — based  on  the  two-material  composite  model. 

According  to  the  formulation  used  here,  the  material  from  which  the  resulting 
structure  is  made  has  optimized  pointwise  properties.  Thus,  the  measure  of  the  objective 
function  corresponds  to  a  lower  bound  of  the  compliance  of  the  design,  relative  to  all 
linearly  elastic  material  properties  orientations,  including  any  possible  configuration  of  a 
composite  material.  In  this  respect,  the  free  elastic  tensor  formulation  produces  a  more 
general  result  than  familiar  methods  founded  on  the  two-material  composite.  On  the 
other,  hand  we  note  that  the  same  general  formulation  may  be  modified  to  provide 
modeling  of  materials  with  manufacturable  microstructure  by  introducing  additional 
constraints  on  the  constitutive  tensor.  The  form  of  such  constraints  should  be  appropriate 
to  reflect  restrictions  on  local  geometry. 

The  external  problem  of  optimal  material  distribution  (P2)  for  given  weighting 
function  is  solved  by  an  iterative  algorithm  based  on  the  optimality  criterion  (16),  taking 
the  first  invariant  of  the  elasticity  tensor  as  the  measure  of  material  properties  density. 
However,  since  the  resulting  design  depends  on  the  choice  of  this  measure,  the  more 
general  problem  is  of  practical  interest.  A  recently  developed  formulation  (Taylor 
(1998)),  where  both  objective  and  cost  are  expressed  in  the  same  basis,  would  be 
convenient  for  this  purpose. 
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Fig.  1 .  Flow  diagram  of  the  computational  procedure 


Fig.2.  Cantilever  beam  subjected  to  transverse  load:  a)  optimal  distribution  of  the  trace  of 
elasticity  tensor  for  continuously  varying  material  properties,  b)-c)  transition  to 
discrete  design  for  slow(b)  and  rapid  (c  )  cut  off  value  adjustment  rate. 


topology  design  (right)  inside 
ter  for  volume  fraction  equal  to: 


Fig.4.  Circular  plate  with  internal  hole  under  uniform  pressure  on  external  edge: 
a)  optimal  solution  for  w  =  1 ,  b)  induced  holes,  c)-g)  transition  to  discrete 
structure,  h)  magnitude  and  directions  of  principal  strains  at  the  final  design. 


Fig. 6.  Final  design  of  3D  cantilever  beam  for  volume  fraction  of  0.2 


Fig/7.  Final  design  of  3D  cantilever  beam  for  volume  fraction  of  0.1 
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Abstract 

The  optimization  of  material  properties  with  the  goal  of  achieving 
minimum  structural  compliance  is  analyzed  in  a  three-dimensional 
setting.  The  formulation  takes  into  account  a  design-dependent  field 
force  with  the  purpose  of  analyzing  problems  where  self- weight  or  cen¬ 
trifugal  forces  have  a  non-negligible  effect.  With  the  introduction  of 
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an  additional  procedure,  optimal  topology  design  is  predicted.  An 
algorithm  based  on  the  optimality  conditions  is  used  to  obtain  nu¬ 
merical  solutions  to  problems  which  illustrate  relevant  effects  of  the 
design-dependent  field  force. 


1  Introduction. 

The  optimal  design  problem  where  the  objective  is  to  achieve  minimum  struc¬ 
tural  compliance  for  a  linearly  elastic  continuum  (or,  equivalently,  maximum 
structural  stiffness)  has  been  used  extensively  as  the  prototypical  example  for 
developing  new  optimization  techniques  and  methods.  Relatively  recently, 
BENDS0E  et.  al  [2]  introduced  a  method  where  the  elasticity  tensor  field  is 
used  as  a  variable  for  the  design  of  a  structure  composed  of  a  general  in- 
homogenous  and  anisotropic  material.  Furthermore,  their  method  predicts 
optimal  material  properties  which  vary  continuously  throughout  the  struc¬ 
ture. 

*********  “black/white”  or  “0/1”  layout  ********** 

Hence,  the  formulation  of  the  design  problem  must  be  modified  in  order 
to  predict  0/1  solutions.  However,  the  modification  should  be  minimal  so 
as  to  retain  the  pratical  relevance  of  the  original  problem.  This  is  precisely 
the  motivation  that  led  Guedes  and  Taylor  [4]  to  introduce  a  piecewise 
constant  weight  function  w  in  an  isoperimetric  constraint  of  the  optimal  de¬ 
sign  problem.  The  idea  was  to  introduce  a  weight  function  taking  only  two 
values,  w  and  w ,  such  that  the  corresponding  solution  to  the  design  problem 
would  also  take  only  two  “values,”  even  though  the  space  of  admissible  de¬ 
sign  variables  is  not  limited  to  such  “values.”  A  stepwise  solution  procedure 
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with  appropriate  adjustement  of  weighting  leads  to  a  result  for  0/1  layout 
design.  The  same  approach  is  considered  here  and  the  main  purpose  is  to 
study  the  effect  of  a  design-dependent  field  force  in  a  three-dimensional  set¬ 
ting.  The  design  of  a  vibrating  elastic  structure  for  maximum  fundamental 
frequency  (see,  e.g.,  Prager  &:  Taylor  [?])  or  the  prediction  of  shape  for 
a  rotating  structure  (Drucker  &;  Shield  [?])  are  examples  of  problems  in 
this  category.  Problems  with  a  field  force  of  this  kind  arise  naturally,  for 
example,  in  the  context  of  the  design  of  structures  where  the  weight  is  taken 
into  account  or  in  the  design  of  elements  belonging  to  a  rotating  structure 
where  a  non-negligible  inertial  load  is  present  (e.g.,  a  turbine,  a  helicopter 
rotor  blade,  etc.).  The  fully  dynamic  problem  is  not  addressed  here;  rather, 
the  inertial  loads  are  introduced ‘as  a  time-independent  field  force  so  that, 
formally,  the  structure  of  an  elastostatic  problem  is  preserved. 

Alternative  methods  to  obtain  a  similar  layout  have  been  proposed  (see, 
e.g.,  the  survey  Rozvany  et.  al  (?]).  One  of  them  is  based  on  the  use  of 
an  artificial  constitutive  power  law  that  penalizes  “intermediate”  values  of 
the  design  variable  (see,  e.g.,  Petersson  and  Sigmund  [7]).  Only  isotropic 
materials  are  used  in  this  model  and  a  (local)  bound  on  the  generalized  deriva¬ 
tive  of  the  design  variable  is  imposed  to  achieve  closure.  Another  technique 
is  to  include  a  term  in  the  objective  functional  that  penalizes  intermediate 
values  of  the  design  variable,  but  independently  of  the  constitutive  law,  by 
using  an  explicit  constraint  (see,  e.g.,  Haber  at  al.  [5]  and  Fernandes  et. 
al  [3]).  This  latter  method  also  includes  a  (global)  bound  on  the  perimeter  of 
the  solid  region  to  achieve  closure.  A  comparative  analysis  of  these  different 
methods  is  beyond  the  scope  of  this  work.  Furthermore,  the  above  mentioned 
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papers  do  not  consider  a  design-dependent  field  force. 

The  basic  optimization  problem  (maximum  structural  stiffness  with  design- 
dependent  field  force)  is  stated  in  Section  2.  The  optimality  conditions  for 
the  design  problem  are  presented  in  Section  3.  Section  4  deals  with  a  version 
of  the  design  problem  where  the  purpose  is  to  generate  a  result  for  0/1  layout 
(i.e.,  topology)  design.  The  numerical  implementation  of  both  the  basic  de¬ 
sign  problem  and  the  0/1  layout  problem  is  developed  in  Section  5  and  some 
numerical  examples  are  presented  in  Section  6.  The  last  section  includes  a 
discussion  and  some  conclusions  regarding  the  effect  of  the  design-dependent 
field  force  in  the  design  of  optimal  structures. 

2  Formulation  of  the  problem. 

2.1  Optimal  design  problem. 

Let  V  be  the  set  of  kinematically  admissible  displacement  fields  and  A  the  set 
of  admissible  design  variables  (fourth-order  symmetric  positive  semi-definite 
tensor  fields  Eijki )•  Consider  an  elastic  body  occupying  a  region  fi,  subject 
to  a  field  force  £><  which  is  assumed  to  depend  on  E^ki-  In  general,  the  field 
force  is  not  necessarily  a  function  of  the  material  properties.  However,  this 
relation  can  be  motivated  in  view  of  the  form  of  the  optimized  structure 
which  is  composed  of  only  two  materials  (i.e.,  0/1  layout).  In  that  case,  the 
elasticity  tensor  field  can  be  viewed  as  a  characteristic  function  that  describe 
the  (disjoint)  domains  occupied  by  each  material,  thus  a  relation  between  bi 
and  E^ki  can  be  established.  In  a  sense,  a  functional  relation  between  the 
field  force  and  the  elasticity  tensor  is  similar  to  an  artificial  constitutive  law 
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where  the  results  are  justified  a  posteriori  only  for  the  optimal  solution. 

Let  W  be  the  work  done  by  the  external  forces,  i.e., 

W(ui,bi(Eijki))  =  [  bi(Eijkt)uidv  +  [  Umda  ,  (1) 

J  n  Jdrit 

where  Uj  is  the  displacement  and  f,-  the  prescribed  surface  traction  (the 
prescribed  displacement  in  dQu  is  zero).  The  potential  energy  is  given  by 
n  =  U  -  W,  where  U  =  5  Jn£ijEijki£kidv  is  the  total  strain  energy.  Recall 
that  the  minimum  value  of  the  potential  energy  is  IT,  =  —  hence,  the 
minimum  structural  compliance  problem  is  stated  as: 

max  min  U(Eijki,Ui)  (Pi) 

Ui£V 

Eijkt€A  . 

0 <£<¥>(£;;«)<?  . 

Jnv{Eijki)dv<R  , 

where  <p(Eijki)  is  a  resource  cost  function,  R  is  an  upper  bound  on  the  total 
resource  and  (p  =  <p(x)  and  Tp  =  Tp{x)  are  given  functions  (independent  of 
Eijki)  satisfying  0  <  <p<Tp.  The  second  explicit  restriction  for  Eijki  is  the 
resource  constraint.  Guedes  and  Taylor  [4]  used  a  weighted  cost  function 
based  on  the  trace  of  E^ki,  i.e.,  ip(Eijki)  =  wEijij,  where  w  is  a  weight 
function.  The  cost  function  considered  here  is  similar,  but  based  on  the 
maximum  eigenvalue  of  Eijki ,  i.e.,  let 

a  =  .max  BijEijkiBki  ,  (2) 

Eij  Bij  —  1 

and  define 


v{Eijki)  =  wa 


(3) 
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where  w  =  w{x)  is  a  positive-valued  weight  function  (an  interpretation  of 
the  cost  function  is  given  in  Section  2.2).  Moreover,  the  functions  <p  and 
Tp  corresponding  to  the  lower  and  upper  bounds  of  the  cost  function  are 
assumed  to  have  the  following  form: 

cp  =  WQ  ,  Ip  =  WQ  , 

where  a  and  a  are  given  positive  constants  satisfying  0  <  a  <  a.  Finally, 
it  is  also  assumed  that  the  field  force  6,  depends  on  a  only,  i.e., 

bi  =  k(Eiju)  =  6i(a)  .  (4) 

It  is  useful  to  consider  discontinuous  weight  functions  (e.g.,  piecewise  con¬ 
stant  functions),  hence  one  has  to  take  into  account  possible  discontinuities 
in  the  design  and  strain  fields  at  points  where  w  is  discontinuous.  However, 
the  displacement  field  is  assumed  to  be  continuous  everywhere. 

2.2  Reformulation  of  the  problem. 

The  design  problem  (Pi)  can  be  reformulated  by  using  a  saddle  point  theorem 
as  proved  in  BENDS0E  et  al.  [2].  This  theorem  can  still  be  applied  with  a 
cost  function  (3)  and  field  force  (4).  In  this  case,  as  shown  in  [8],  the  cost 
function  (3)  provides  a  positive  definite  optimal  material  of  the  form 

E°jU  =  (a  -  a)-^-  +  Ifl  +  M>t)  .  (5) 


The  corresponding  stress  field  is 
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and  the  potential  energy  is 

n0(a,  Ui)  =  /  ^oc£ij£ijdv  -  W0{ui,bi{cc))  .  (7) 

J  ci' 

The  optimal  material  characterized  by  (5)  is  an  anisotropic  (orthotropic) 
material.  However,  due  to  the  coupling  between  the  design  and  elastostatic 
problems  and  in  view  of  (6),  the  optimal  material,  restricted  to  the  opti¬ 
mal  strain,  behaves  like  an  isotropic  material  with  zero  Poisson’s  ratio  and 
Young  modulus  equal  to  a.  The  optimal  material  can  be  decomposed  into  an 
isotropic  part  characterized  by  a  Young’s  modulus  E  =  a  and  zero  Poisson’s 
ratio  and  an  anisotropic  “reinforcement”  aligned  with  the  optimal  strain  field 
c ij  and  characterized  by  a  unique  non-zero  value  a  —  a  (i.e.,  essentially  it 
depends  on  the  maximum  eigenvalue  a  of  Eijki)-  The  isotropic  part  of  Efjkl 
is  fixed  and  the  problem  reduces  to  identify  the  anisotropic  reinforcement  in 
order  to  optimize  the  structure.  The  cost  function  measures  the  (weighted) 
amount  of  reinforcement  required  for  the  optimal  structure.  Thus,  the  opti¬ 
mal  design  problem  can  be  reformulated  as 

max  min  n0(a,Ui)  .  (P2) 

a 

0<Q  <Q<Q  , 
ti/adv<R  , 

At  points  where  a  =  a,  the  interpretation  is  that  no  reinforcement  is  re¬ 
quired. 
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3  Optimality  conditions. 


Consider  the  Lagrangian  L  given  by 


L(a,  Ui ;  A,  A,  A)  =  n0(ar,  uj  +  J  [A(q  -  a)  -  A(a  -  a )]  dz; 


4-  A  jy  wadv  -  R  j  , 


where  A  =  A(x),A  =  A(x)  and  A  are  Lagrange  multipliers.  The  presence 
of  a  (possibly)  discontinuous  weight  function  needs  to  be  taken  into  account 
since  the  design  variable  or  and  the  strain  field  would  also  be  discontinuous  at 
the  same  points.  Let  S  be  the  set  of  such  points  of  discontinuity.  Arbitrary 
variations  of  L  with  respect  to  U{  provide  the  equilibrium  equations.  Since 
U{  is  continuous,  then  so  is  5ui\  in  particular,  =  0  on  the  surface  of 
discontinuity  S,  where,  for  any  function  /,  [/]  =  f  +  —  f  ,  f±  being  the 
values  of  j  on  each  side  of  S.  For  given  or,  the  local  form  of  the  elastostatic 
problem  for  the  optimal  strain  at  points  where  a  and  dj  are  smooth  is 

(a£ij)j  +  bi(a)  =  0  infi  — S, 

OijTij  =  ctSijTij  =  t{  on  dftt  ,  ?  (9) 

U{  =  0  on  dQu  , 

/ 

where  rij  is  the  outward  unit  normal  vector  to  the  boundary  dFl .  At  points 
where  a  and  are  discontinuous,  one  has 


J [aeijSu^Uida  =  J  {[ar^j'J^^i)  +  (ors.j)^^]}  da  =  0 


where  Uj  is  the  normal  vector  pointing  towards  the  “+”  side  of  the  surface  of 
discontinuity  S  and,  for  any  function  /,  (/)  =  i(/+  +  /“).  Moreover,  since 
the  integration  by  parts  does  not  involve  <56,-,  possible  discontinuities  of  6;  do 
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not  affect  the  result.  Therefore,  since  [<fu;J  =  0,  the  local  form  of  the  above 
equation  is 

[o’Sijji/j  =  0  on  S  .  (10) 

(This  condition  represents  the  continuity  of  the  traction.)  Let  V  be  the 
gradient  of  L  with  respect  to  a,  i.e.,  5L[5a\  =  fn  L'Sadv  and  let  b\  be  the 
gradient  of  6,-  with  respect  to  a;  the  gradient  V  can  be  expressed,  Vx  £  ft  — S, 
by 

L'  =  b^Ui  —  \tijZ\j  —  A  +  A  +  Au;  .  (11) 

Therefore,  the  Karush-Kuhn-Tucker  optimality  conditions  are,  Vx  £  ft  —  S, 

L‘  =  0  , 

A(a  —  a)  =  0  , 

A(a  —  a)  =  0  , 

A  {/n  wadv  —  i?}  =  0  , 

At  points  on  S,  the  jump  condition  [L‘ J  =  0  is  trivially  satisfied.  For  a  more 
specific  formulation  of  the  optimality  conditions  one  can  use  the  following 
domains: 

D  =  {x  £  Cl  |  a(x)  =  a  }  , 
ft  =  {x  £  Q  |  a(x)  =  a  }  , 
fti  =  (x  £  ft  |  a  <  a(x)  <  a  }  . 


A  >  0  , 
A  >  0  , 
A  >  0  . 


(12) 


(13) 
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Hence,  since  A  =  0  for  x  €  Ft,  A  =  0  for  x  6  Q  and  A  =  A  =  0  for  x  €  Qt-,  the 
optimality  conditions  (12)  become 


Aw  —  2 Eij^ij  , 

x  6  Cli 

Aw  X  —  <2,Eij£ij  , 

A  >  0  , 

i6fi, 

Aw  +  A  =  \eijEij  -  b'jUi  , 

A  >  0  , 

x  €  Cl  , 

A  (fn  wadv  -  Rj  =  0  , 

A  >  0  , 

x  €  Cl  . 

In  Section  4,  a  class  of  discontinuous  functions  w  will  be  introduced  such  that 
the  set  of  points  of  discontinuity  S  can  lie  either  inside  or  on  the  boundary 
of  Di  or  on  the  boundaries  between  Q  and  Cl.  In  any  case,  for  that  class  of 
functions,  the  jump  condition  is  trivially  satisfied,  i.e.,  from  (14), 

AM  =  ~  b'iUi  +  A  -  A]  on  S  . 

If  w  is  continuous,  then  the  jump  conditions  do  not  apply.  Moreover,  if 
b'i  =  0  then  one  can  prove  that  A  >  0  (and  hence  that  the  resource  constraint 
is  satisfied  as  an  equality,  see  [4]).  However,  this  is  not  necessarily  true  if 
b\  ^  0  (i.e.,  it  is  possible  that  A  =  0  and  the  problem  (P2)  has  a  non-trivial 
solution).  Before  closing  this  section,  it  is  interesting  to  observe  that,  from 
(14),  the  term  b'^i  +  Aw  is  necessarily  non-negative  in  f2,Uf2  (even  if  A  =  0), 
but  it  could  have  any  sign  in  Cl.  Thus,  if  b\ui  +  Aw  <  0,  then,  necessarily, 
this  can  only  happen  at  a  point  where  the  material  is  at  its  upper  bound. 
This  fact  will  be  used  in  Section  5. 
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4  Optimal  design  problem  with  discrete  val¬ 
ues. 

As  mentioned  in  the  introduction,  as  an  additional  design  requirement  it  is 
useful  to  identify  solutions  to  the  problem  (P2)  where  the  design  variable  is 
restricted  to  only  two  values:  a  and  a .  The  resulting  structure  is  charac¬ 
terized  by  elasticity  tensors  {Sik^ji  +  $u$jk)  and  (a  —  a  )Eij£ki/{smnemn)  + 
ko.  {fiikfiji  +  fiufijk)  respectively.  The  region  occupied  by  the  first  material 
does  not  have  “reinforcement,”  whereas  the  complementary  region  has  a 
(reinforced)  material  with  fixed  magnitude  a  and  variable  orientation  (the 
magnitude  of  E,jki  is  measured  here  in  terms  of  its  maximum  eigenvalue;  its 
Frobenius  norm  is  ^/EijkiEijki  =  \/a (a  —  a)  +  6a2).  For  the  specific  case 
where  a  and  a  are  chosen  such  that  a  a,  then  the  two  materials  are 
interpreted  as  “void”  and  “full  material.”  In  general,  this  corresponds  to 
a  suboptimal  solution  compared  to  one  where  o  varies  continuously  (which 
can  be  obtained,  e.g.,  with  a  weight  function  w  =  1).  However,  it  is  an 
optimal  solution  for  some  other  weight  function  and  satisfies  the  additional 
requirement  of  taking  only  two  values. 

The  use  of  a  weight  function  is  intimately  related  to  the  resource  con¬ 
straint.  In  general,  using  only  local  bounds,  one  cannot  rule  out  trivial  so¬ 
lutions.  The  resource  constraint  is  used  to  optimally  distribute  the  material 
properties.  However,  the  distribution  of  material  properties  that  satisfy  the 
resource  constraint  is  accomplished  via  a  single  (constant)  parameter,  i.e., 
via  the  Lagrange  multiplier  A.  In  a  sense,  the  use  of  a  weight  function  w  is 
an  extension  of  this  idea,  except  that  w  acts  locally  via  the  product  Aw  as 
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shown  by  (14).  Hence,  the  effect  of  the  Lagrange  multiplier  can  be  modified 
locally  by  a  proper  use  of  a  weight  function. 

Let  F  be  the  class  of  piecewise  constant  weight  functions  defined  as 


w  G  F  w  = 


w  in  Cl™  , 
w  in  Cl™  , 


(15) 


where  the  domains  Cl™  and  Cl™  (to  be  defined  later)  satisfy  Cl™  U  Cl™  =  Cl 
and  Cl™  fl  =  0.  The  superscript  reflects  the  fact  that  different  domains 
correspond  to  different  functions  w.  Suppose  that  for  a  given  weight  function 
the  corresponding  solution  to  (P2)  is  a  (in  general,  not  a  piecewise  constant 
function).  Recall  that  the  domain  Qit  as  defined  by  (13),  represents  the 
points  where  the  design  variable. is  between  the  lower  and  upper  bounds  a 
and  a .  As  an  additional  objective  one  looks  for  a  function  w.  €  F  in  the 
following  problem: 


inf  in,  I  ,  '  (P3) 

t v£F 

where  Cl,  is  defined  as  an  implicit  function  of  w  by  (13)3,  (15)  and  the  cor¬ 
responding  solution  a  to  (P2).  The  functional  dependence  of  fit-  on  w  is 
somewhat  complex  but  numerical  experiments  have  shown  that  there  are 
(numerical)  approximations  of  a  function  u/„  such  that  |Q,-|  «  0.  The  result¬ 
ing  structure  satisfies  Cl  =  Cl  U  Cl  and  Q  fl  Cl  =  0.  However,  the  uniqueness, 
and,  in  fact,  the  existence  of  w,  in  the  general  case,  have  not  been  established. 
Nevertheless,  the  purpose  here  is  to  find  at  least  one  numerical  approxima¬ 
tion  to  such  a  function  w..  The  domains  Cl™  and  Cl™  of  (15)  are  now  defined 
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as 

fi!  =  {i  6  fi  I  a(i)  <  ac}  , 

(16) 

=  {z  €  ft  |  a(x)  >  ac}  , 

where  the  constant  ac,  a  <  ac  <  a ,  is  a  prescribed  “cutoff’  value.  The  in¬ 
terpretation  of  the  weight  function  w  follows  from  this  definition:  the  method 
assigns  a  weight  w  to  regions  of  “weak”  reinforcement  (a  <  ac)  and  a  weight 
ui  otherwise. 

In  order  to  solve  the  problem  (P3),  an  iterative  procedure  is  used  since,  a 
priori,  the  domains  and  cannot  be  defined  without  knowledge  of  the 
optimal  a  which,  in  turn,  depends  on  w.  The  details  of  the  algorithm  will  be 
given  in  Section  5,  but  it  is  interesting  to  point  out  a  few  facts  here:  suppose 
that  there  exists  a  minimizer  w,  of  (P3)  and  assume  that  min,^/?  |fij|  =  0, 
then 

hence  the  region  of  weight  w,  coincides  with  the  region  without  rein¬ 
forcement  and  the  region  of  weight  w ,  coincides  with  the  reinforcement. 
From  this  point  of  view,  the  function  w.  is  related  to  the  characteristic  func¬ 
tions  of  £2  and  Cl.  Moreover,  recall  that  the  resource  constraint  is  viewed  as 
a  global  bound  on  the  “cost”  of  reinforcement.  If  the  resource  constraint  is 
active  for  the  weight  function  w.  (i.e.,  A  >  0),  one  has,  from  (14)4, 

aw  |Q|  +  aw  |ft|  =  R  , 

thus,  for  an  appropriate  interpretation  of  the  regions  with  and  without  rein¬ 
forcement,  one  should  choose  w  and  w  such  that 
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If  Q.  |0I  /a  | ft  |  is  bounded  above  and  below,  then  a  choice  such  as  w  «  w 
seems  appropriate.  However,  based  on  the  optimality  condition  (14)o,  this 
choice  is  not  consistent  with  a  =  a  in  ft.  This  can  be  seen  by  assuming 
b'i  =  0,  hence  the  multiplier  A  would  not  be  positive  in  ft  if  w  is  “small”. 
Nonetheless,  suppose  that  a  «  a  (which  corresponds  to  the  void/full  ma¬ 
terial  case)  and  assume  that  |ft|  /  |ft|  is  bounded  above  and  below,  then 
a,a,w  and  w  should  satisfy 

q  «a  ,  w  w  ,  aw  <  aw  , 

e.g.,  ifa/5  =  10-P,  w  /w  =  10-?,  aw /aw  =  then  one  should  take 

p  >  0,q  >  0  and  p  -  q  >  0  sufficiently  large,  which  gives  R  ~  a|ft|  (with 
w  =  1).  The  above  choice  does  not  contradict  the  optimality  conditions, 
therefore,  the  use  of  a  weight  function  will  be  limited  to  the  void/full  material 
case.  The  resource  cost  function,  which  is  the  product  of  the  weight  function 
and  the  design  variable  a,  would  be  small  when  a  =  d  and  large  when 
a  =  a .  Hw  and  w  are  of  the  same  order  of  magnitude,  then,  in  general,  the 
corresponding  structure  would  not  be  void/full  material.  The  method  does 
not  require  to  set  w  =  1,  but  this  choice  provides  the  same  interpretation  of 
the  resource  constraint  as  in  the  case  when  no  weight  function  is  used. 

5  Numerical  implementation. 

5.1  Overview. 

The  algorithm  has  two  main  elements:  the  first  one  deals  with  the  basic  op¬ 
timization  problem  (P2)  for  a  given  (fixed)  weight  function  and  the  second 
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one  corresponds  to  a  sequence  of  basic  optimization  problems  solved  for  dif¬ 
ferent  weight  functions  with  the  purpose  of  identifying  a  minimizer  w.  in  the 
problem  (P3). 

In  order  to  obtain  a  numerical  solution  of  the  optimal  design  problem 
(P2),  a  method  based  on  the  optimality  conditions  (14)  is  used.  One  can 
compute  a  numerical  sequence  {an}  of  the  design  variable  such  that  each  an 
is  chosen  within  the  constraint  space.  This  is  enforced  via  the  local  and  global 
Lagrange  multipliers.  In  order  to  minimize  W0  in  (P2),  the  corresponding 
displacement  field  Uj  (for  a  given  a)  has  to  be  obtained  as  the  solution  of  (9). 
This  can  be  achieved  using  a  standard  finite  element  code  (in  this  case  the 
program  ANSYS  was  used). 

For  the  problem  (P3),  the  basic  optimization  problem  (P2)  is  solved  with 
a  given  weight  function  and  a  new  weight  function  is  defined  from  (15)  and 
(16)  using  the  new  distribution  of  material  properties  a.  A  new  optimization 
problem  (P2)  is  solved  with  the  updated  weight  function  until  a  0/1  layout  is 
identified.  Section  5.2  deals  with  the  basic  optimization  problem  (P2)  while 
the  0/1  procedure  for  problem  (P3)  is  described  in  Section  5.3. 

5.2  Basic  optimization  problem. 

The  details  of  the  algorithm  are  as  follows  (a  concise  version  is  given  af¬ 
terwards).  Suppose  that  a  given  approximation  an  has  been  identified.  To 
obtain  an+i,  first  compute  u"  (and  hence  £?•)  by  solving  problem  (9)  with 
a  =  an.  To  remain  in  the  constraint  space  at  step  n  + 1  while  minimizing  the 
objective  functional,  one  needs  to  iterate  in  order  to  satisfy  all  constraints 
simultaneously.  (The  superindex  k  will  be  used  in  connection  with  this  iter- 
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ative  procedure.)  The  method  proposed  here  is  similar  to  the  one  used  by 
Bendsoe  et  al.  [1]  except  that  it  is  not  assumed  that  the  global  Lagrange 
multiplier  is  positive  since  the  resource  constraint  might  become  active  from 
step  n  to  step  n  +  1  or  might  not  be  active  at  all  (this  is  due  to  the  presence 
of  a  design-dependent  field  force).  First,  update  an  based  on  the  optimality 
conditions  without  including  the  local  constraints,  i.e.,  from  (14),  define 


an+i  = 


1 

2_y  v 


IP 


max{7,  +  Aj+,t»} 


(17) 


where  7  >  0  is  a  “small”  number  and  p  is  an  appropriately  chosen  value 
to  improve  convergence  (usually  0  <  p  <  1).  For  k  =  0,  one  can  take 
A®+1  =  An  which  corresponds  to  the  last  value  of  A  at  step  n.  The  form  of 
the  denominator  follows  from  the  remark  at  the  end  of  Section  3:  if  (6”)'u"  + 
An-i-i^  <  0,  then,  setting  it  equal  to  a  small  number  will  force  it  to  its  proper 
value,  i.e.,  the  upper  bound.  The  power  p  is  used  to  avoid  large  fluctuations  of 
the  design  variable.  To  complement  the  effect  of  p,  and  additional  restriction 
can  be  implemented  by  imposing  a  “move  limit”  £  as  follows: 


*n+l  =  < 


O' 


Wi 


(1  -  C)an 

^  <  (1  “  C)“n  , 

6n+l 

if  (1  -  Qan  <  d*+1  <  (1  +  C)an  . 

(18) 

(i  +  CK 

if  o£+1  >  (1  +  CK  • 

e  local  constraints,  i.e.,  for  k  >  0, 

QL 

if  A*+1  <  a  , 

<*n+l 

if  QL  <  <5n+i  <  5  . 

(19) 

a 

if  Q*+i  >  5  • 

Observe  that,  without  a  denominator  as  the  one  used  in  (17),  the  value  of 
when  (6”)'u”  +  A*+1u/  <  0  would  be  a  instead  of  its  correct  value  a. 
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If  0  <  (b’i)'Ui  +  Rn+iw  <  7>  then  the  denominator  in  (17)  has  the  same 
effect,  however,  this  does  not  preclude  q*+1  to  take  values  different  than  a 
since  the  numerator  could  be  small  as  well.  In  order  to  satisfy  the  resource 
constraint  (either  as  an  equality  or  an  inequality),  a  penalty  method  is  used. 
The  Lagrangian  is  augmented  by  adding  the  term 

577  wadv  -  F^j  , 

where  77  >  0  is  the  penalty.  Strictly  speaking,  the  optimality  conditions  (14) 
should  be  modified  to  account  for  this  term,  however  this  is  taken  care  of  in 
the  numerical  procedure.  Define 

74+1  =  [  ,  .  A*  =  -  R  ,  (20) 

J  n 

and  update  the  Lagrange  multiplier  as  follows: 

An+i  =  max  {0,  A*+1  +  77A*}.  (21) 

Observe  that  this  includes  the  case  An+i  =  0  if  the  resource  constraint  is 
not  active.  With  this  new  multiplier,  one  can  compute  from  (17),  then 
from  (18)  and  impose  the  local  constraints  (19)  to  obtain  as  before. 
The  process  converges  when  |A^  —  A£+1|  <C  1,  in  which  case  1 

if  the  constraint  is  active.  At  the  end  of  this  iterative  procedure  one  should 
have  an  appropriate  approximation  an+i  which  satisfies  all  constraints  based 
on  the  previous  step’s  strain  field  £?•.  The  procedure  to  obtain  a  numerical 
solution  of  (P2)  can  be  summarized  in  the  form  of  the  following  algorithm: 

A.  Initialization. 


1:  Set  q0  =  ce  (Given  initial  guess). 
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2:  Obtain  u°,  £?•  from  (9). 

3:  Set  Ao  =  A  (Given  initial  guess;  could  be  zero). 

B.  Update:  orn  — >  an+i- 

1:  Set  A°+1  =  A„  and  compute  a°+1,  q°+1  and  q°+1  from  (17),  (18) 
and  (19). 

2:  For  k  >  0:  compute  R^+i  and  A*  from  (20)  and  update  A£+1  -4 
A*tj  from  (21)  and  then  o£+1  -4  a**}  from  (17),  (18)  and  (19). 

3:  Check  if  |A^{  -  A*+1|  <  1,  in  which  case  take  an+i  =  a*i{, 
otherwise  repeat  step  3  until  convergence. 

C.  Convergence  (objective  functional):  check  if  |WJl+1  —  kF0n|/|kkr0n|  -C  1, 
otherwise  repeat  B  until  convergence.  Alternatively,  one  could  check 
the  optimality  conditions  (13)  in  some  average  sense  (e.g.,  L2-norm). 

5.3  Discrete  (0/1)  layout. 

The  algorithm  presented  here  is  related  to  the  problem  (P3)  as  described 
in  Section  4.  Recall  that  in  this  case  the  values  of  o,q,w  and  w  cannot 
be  chosen  arbitrarily.  In  order  to  obtain  a  discrete  layout,  a  sequence  of 
optimization  problems  (P2)  are  solved.  For  the  first  optimization  problem, 
a  uniform  weight  function  is  prescribed:  the  initial  cutoff  value  is  o°c  =  0, 
hence,  from  (16), 

wq  —  w  Vx  6  n . 

At  step  /,  solving  the  optimization  problem  with  a  weight  function  wi,  using 
the  procedure  described  in  Section  5.2,  provides  the  corresponding  optimal 
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design  a*.  The  weight  function  at  step  /  + 1  is  obtained  from  step  l  as  follows: 
the  cutoff  value  is  incremented  by  a  (possibly  constant)  value  Actc, 

=  o£  +  Aac  , 

and  new  domains  and  +1  are  defined  from  (16)  using  a-j  (hence  wi+i 
is  determined  from  (15)).  For  practical  purposes,  one  can  choose  a  sufficiently 
large  number  of  iteration  steps,  say  nj,  and  define  a  constant  increment 

a  —  a 

&ac  = - , 

Til 

although  a  variable  adjustment  is  sometimes  more  efficient.  The  process 
converges  if  |nj+1|/|fl|  <  1.  A  proof  of  convergence  is  not  available  but 
numerical  experiments  show  that  it  is  possible  to  attain  convergence  with 
Acte  small  enough.  If  Og+1  >  a  but  |fi{+l|  0,  then  the  method  fails.  The 

procedure  is  summarized  in  the  following  algorithm: 

t 

A.  Initialization. 

1:  Set  q°  =  0,wo  =  w . 

2:  Solve  optimization  problem  (P2)  as  described  in  Section  5.2  with, 
e.g.,  d  =  constant  (initial  guess  for  a),  and  obtain  the  correspond¬ 
ing  solution  oro  (the  subindex  in  o;0  refers  to  the  solution  of  (P2) 
with  w  =  wo,  not  to  the  initial  guess). 

B.  Update:  wi  — >  wt+\. 

1:  Set  a(.+1  =  a[  +  Aac. 

2:  Obtain  wi+i  from  (15)  and  (16)  using  Q;  and  a(.+1. 
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3:  Solve  optimization  problem  (P2)  with  a  =  and  obtain  at+ 1. 

C.  Convergence  (0/1  layout):  check  if  |fi[+1|/|f5|  •C  1,  otherwise  repeat  B 
until  convergence. 

The  application  of  this  algorithm  to  an  example  problem  is  described  in  the 
next  section. 


6  Example  applications:  the  optimal  rotating 
structure. 

As  an  illustrative  three-dimensional  example,  consider  the  case  of  the  (stiffest) 
design  of  a  rotating  structure  occupying  a  design  domain  f 2  and  rotating  at 
a  constant  angular  velocity  u.  Let  {e-,c?,c?}  be  a  cylindrical  orthonormal 
basis  where  e?  is  aligned  with  the  axis  of  rotation  and  e},ef  are  fixed  with 
respect  to  fi.  The  domain  Q  is  taken  as  a  prismatic  bar  of  unit  lenght  and 
rectangular  cross  section  of  area  h 2.  Suppose  that  the  mass  density  p  is 
related  to  the  design  variable  a  linearly,  i.e., 

p  =  Ka  , 

where  k  is  a  positive  constants.  In  this  case,  the  inertial  force  can  be  modeled 
as  the  body  force 

bi(a)  =  u7K.rae]  ,  (22) 

where  r  is  the  orthogonal  distance  from  a  given  point  in  Q  to  the  axis  of 
rotation.  The  side  attached  to  the  axis  of  rotation  is  modeled  as  a  clamped 
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end.  The  remaining  sides  of  the  structure  are  subjected  to  the  following 
loads: 


Case  1:  shearing  load: 

{Toef  for  r  =  l  (Uniform  shear  stress) 
0  otherwise 


Case  2:  torsional  load: 

/ 

-r0ef  for  r  =  1  , 
U  =  \  T0e?  for  r  —  1  , 

0  otherwise 


x  e  [-h,  -h  +  e]  , 
x  e  [h  -  e,  h]  , 


For  the  loading  Case  2,  e  is  a  small  number,  so  that  two  opposite  uniform 
shear  loads  are  applied  on  strips  pf  width  e.  Recall  that  h  is  the  width  of  the 
square  cross-section  and  x  is  measured  in  the  e?  direction. 

For  illustration  purposes,  the  corresponding  optimal  distribution  of  ma¬ 
terial  properties  (where  a  =  a)  is  shown,  for  Case  1,  in  Figures  1,  2  and  3 
with  and  without  a  field  force  (i.e.,  with  and  without  rotation).  The  design 
without  a  field  force  is  symmetric  with  respect  to  the  mid-plane  ej ,  ef  but 
by  taking  into  account  the  field  force  the  structure  is  reinforced  on  the  top 
(which  is  in  compression)  and  is  weakened  on  the  bottom  (which  is  in  ten¬ 
sion).  Figures  4  and  5  correspond  to  the  design  for  loading  Case  2.  Observe 
that  the  interior  of  the  structure  is  essentially  empty;  in  the  case  without  a 
field  force,  it  reproduces  a  well-known  optimal  design:  a  prismatic  body  with 
circular  cross-section  (within  the  limitations  of  a  coarse  mesh).  However,  the 
design  with  a  field  force  is  similar  but  it  is  not  prismatic:  it  is  tappered; 
thicker  close  to  the  axis  of  rotation  and  thinner  towards  the  end  where  the 
load  is  applied. 
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7  Discussion  and  conclusions. 

The  differences  between  the  optimal  designs  with  and  without  a  field  force, 
as  shown  in  the  examples  of  Section  6,  can  be  traced  back  to  the  same  term 
in  the  optimality  conditions  (14):  if  b\ui  <  0,  the  design  is  locally  stronger 
(compared  to  the  case  without  a  field  force,  i.e.,  bi  =  0).  The  converse  is 
true  when  >  0.  In  Case  1  (shearing  load),  the  top  (resp.  bottom)  of  the 
structure  is  in  compression  (resp.  tension),  hence  b\ui  <  0  (resp.  >  0) 
and  this  accounts  for  the  features  in  the  design.  In  Case  2  (torsional  load)  the 
situation  is  slightly  different  than  in  Case  1  since  >  0  everywhere,  yet, 
compared  to  the  case  without  field  force,  some  sections  are  weaker  but  other 
are  stronger.  In  this  case,  the  term  is  greater  towards  the  end  of  the 
structure  where  the  surface  loads  are  applied.  Since  in  this  case  the  (global) 
resource  constraint  is  active,  there  is  a  fixed  amount  of  material  that  is  used, 
hence  the  structure  becomes  thicker  closer  to  the  axis  of  rotation  where  6'u,- 
is  smaller,  as  illustrated  in  Figure  5. 

Although  some  technical  issues  remain  open,  the  method  shown  here 
provides  a  practical  way  of  identifying  a  0/1  layout  with  an  account  of  a 
design-dependent  field  force.  Furthermore,  the  method  is  relatively  easy 
to  implement  for  three-dimensional  problems  since  standard  finite  element 
programs  can  be  used  for  this  purpose.  As  a  closing  remark,  it  is  worth 
noting  that  the  use  of  a  weight  function  is  intimately  related  to  the  resource 
constraint  since,  if  this  constraint  is  not  active  (A  =  0),  the  weight  function 
does  not  affect  the  solution  of  the  optimal  design  problem.  Hence,  in  that 
case,  one  cannot  expect  to  obtain  a  0/1  layout  using  this  method. 


February  11,  1998  15:52  DRAFT 


23 


References 

[1]  Bendsoe,  M.,  Diaz,  A.  and  Kikuchi,  N.:  Topology  and  gener¬ 
alized  layout  optimization  of  elastic  structures,  in:  Topology  Design 
of  Structures,  M.  Bendspe  and  C.  A.  Mota  Soares  (eds.),  Kluwer, 
NATO  ASI  series,  227,  pp  159-205,  1993. 

[2]  Bendsoe,  M.  P.,  Guedes  J.  M.,  Haber  R.  B.,  Pedersen 
P.  and  Taylor,  J.  E.:  An  analytical  model  to  predict  optimal 
properties  in  the  context  of  optimal  structural  design,  Journal  of 
Applied  Mechanics,  61,  pp  930-937,  1994. 

[3]  Fernandes,  P.  R.,  Guedes,  J.  M.  and  Rodrigues,  H.,  Topol¬ 
ogy  optimisation  of  3D  linear  elastic  structures,  Publication  pending. 

[4]  Guedes,  J.  M.  and  Taylor,  J.  E.,  xxx,  xxx,  xx,  pp  xxx-yyy, 
19zz. 

[5]  Haber,  R.  B.,  Jog,  C.  S.  and  Bendsoe,  M.  P.:  A  new  approach 
to  variable-topology  shape  design  using  a  constraint  on  perimeter, 
Structural  Optimization,  11,  pp  1-12,  1996 

[6]  Jog,  C.  S.,  Haber,  R.  B.  and  Bendsoe,  M.  P.:  A  displacement- 
based  topology  design  method  with  self-adaptive  layered  materials, 
in:  Topology  Design  of  Structures,  M.  Bendspe  and  C.  A.  Mota 
Soares  (eds.),  Kluwer,  NATO  ASI  series,  227,  pp  219-238,  1993. 


February  11,  1998  15:52  DFLAFT 


24 


[7]  Petersson,  J.  and  Sigmund,  0.:  Slope  constrained  topology  opti¬ 
mization,  To  appear  in  International  Journal  of  Numerical  Methods 
in  Engineering. 

[8]  Turteltaub,  S.  and  Washabaugh,  P.,  xxx,  xxx,  Publication 
pending. 


February  11,  1998  15:52  DRAFT 


25 


Figure  1:  Case  1  (shear  force):  isometric  plot  of  the  optimal  design  without 
a  field  force  (top)  and  with  a  field  force  (bottom).  For  clarity,  only  half  of 
the  design  is  shown  (i.e.,  for  x  €  [—  h,  0]). 


No  rotation 


Figure  2:  Case  1  (shear  force):  right  view  of  the  optimal  design  without  a 
field  force  (top)  and  with  a  field  force  (bottom). 
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No  rotation 


Rotation 


Figure  3:  Case  1  (shear  force):  top  view  of  the  optimal  design  without  a  field 
force  (left)  and  with  a  field  force  (right). 
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Figure  4:  Case  2  (torsional  force):  isometric  view  of  the  optimal  design 
without  a  field  force  (top)  and  with  a  field  force  (bottom).  For  clarity,  only 
half  of  the  design  is  shown  (i.e.,  for  x  €  [—  h,0]). 
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CO 

Torsional  load  I 


r  -  0.2 


No  rotation 


Rotation 


Figure  5:  Case  2  (torsional  force):  cross-sectional  view  the  of  optimal  design 
(for  t  —  0.2)  without  a  field  force  (left)  and  with  a  field  force  (right). 
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Abstract  A  method  is  presented  for  the  prediction  of  optimal  configurations  for  two-material 
composite  continuum  structures.  In  the  model  for  this  method,  both  local  properties  and  topology 
for  the  stiffer  of  the  two  materials  are  to  be  predicted.  The  properties  of  the  second,  less  stiff 
material  are  specified  and  remain  fixed.  At  the  start  of  the  procedure  for  computational  solution, 
material  composition  of  the  structure  is  represented  as  a  pure  mixture  of  the  two  materials.  This 
design  becomes  modified  in  subsequent  steps  into  a  form  comprised  of  a  skeleton  of  concentrated 
stiffer  material,  together  with  a  nonoverlapping  distribution  of  the  second  material  to  fill  the 
original  domain.  Computational  solutions  are  presented  for  two  example  design  problems.  A 
comparison  among  solutions  for  different  ratios  of  stiffness  between  the  two  materials  gives  an 
indication  of  how  the  distribution  of  concentrated  stiffer  material  varies  with  this  factor.  An 


example  is  presented  as  well  to  show  how  the  method  can  be  used  to  predict  an  efficient  layout  for 
rib-reinforcement  of  a  stamped  sheet  metal  panel. 


1.  Introduction 

We  address  the  problem  of  how  to  predict  the  optimal  configuration  of  a  continuum  structure 
composed  of  two  distinct,  linearly  elastic  materials.  The  problem  is  cast  in  a  form  that  has  the 
stiffer  of  the  component  materials  treated  as  ‘designate’,  while  properties  of  the  second,  less-stiff 
material  are  taken  to  be  specified  and  fixed.  The  goal  in  the  treatment  of  this  design  problem  is 
twofold:  i)  to  determine  the  local  properties  of  the  concentrated  stiffer  material  and  ii)  to 
determine  its  the  optimal  topology,  imbedded  in  the  second  material  which  occupies  all  the 
remainder  of  the  original  domain  of  the  structure.  The  intention  is  that  this  result  simulates  an 
optimally  reinforced  two-material  composite  continuum  structure.  An  interpretation  of  the  model 
for  computational  solution  is  applied  to  produce  example  results  in  2D  and  3D  showing  the" form 
of  such  designs.  As  another  type  of  application,  the  modeling  technique  is  applied  to  design  the 
optimal  layout  for  rib-reinforcement  in  a  stamped  sheet  metal  panel. 

The  analytical  formulation  for  the  composite  design  problem  is  based  on  the  concept  that  has  a 
designable  continuum  material  represented  (for  linearly  elastic  material)  by  its  constitutive  tensor. 
The  paper  by  Bendsoe  et  al.  (1994)  evidently  is  the  first  example  where  the  arbitrary  tensor- valued 
function  representing  material  properties  is  treated  directly  as  the  design  variable.  [A  quite 
different  approach  to  represent  variable  material  property  within  a  design  problem  is  described  in 
Jacobs  et  al.  (1997)].  In  that  presentation  as  well  as  subsequent  other  applications  (see  e.g. 
Bendsde  et  al.  (1995,  96);  Bendspe  (1995)],  the  design  objective  was  to  minimize  compliance  and 
the  global  cost  (isoperimetric)  constraint  was  expressed  in  terms  of  the  trace  or  the  second  order 
invariant  of  the  modulus  tensor.  The  present  formulation  is  stated  for  the  same  objective,  and  for 
the  cost  constraint  based  on  the  trace  measure  of  the  modulus  tensor. 

Results  from  the  cited  earlier  design  formulations  are  comprised  of  a  prediction  of  the  local 
form  of  the  optimal  material  (  a  zero-shear  stiffness,  orthotropic  material  in  the  case  of  single¬ 
purpose  design  ),  together  with  the  (continuously  varying)  distribution  of  the  trace  of  the  modulus 
tensor,  the  latter  representing  a  measure  of  merit  of  the  material.  In  the  present  setting,  the 
corresponding  result  has  the  form  of  an  optimal  distribution  of  a  pure  mixture  of  the  two 


constituent  materials,  having  effective  modulus  equal  to  the  sum  of  the  constituent  moduli.  In 
order  to  be  able  to  determine  the  design  for  a  physically  realistic  composite  material,  i.e.  a 
structural  composite  where  the  two  materials  are  distinct,  a  recently  described  technique  [Guedes 
and  Taylor  (1997a,b),  Pawlicki  et  al.  (1998)  for  example  applications  in  3D]  for  the  computational 
prediction  of  optimal  topology  is  applied.  The  technique,  which  differs  sharply  from  the  familiar 
approach  [see  e.g.  Bends0e  (1995),  Olhoff  et  al.  [1997],  amounts  to  a  step  by  step  computational 
procedure,  making  use  of  the  above  described  continuously  varying  design  as  a  starting  point.  A 
finite  sequence  of  repeated  solutions  to  the  original  design  problem,  each  with  stepwise,  ordered 
modification  to  the  ‘unit  cost  distribution’,  leads  to  the  final  result  in  the  form  of  a  composite 
having  the  two  materials  appear  as  effectively  distinct  but  mechanically  combined.  (A  distinctly 
different  example  where  ’overlap  of  materials’  is  addressed  is  reported  in  Rozvany  et  al.  (1982).) 

To  summarize  the  contents  of  the  paper,  the  model  for  prediction  of  the  optimal  composite 
continuum  structure  is  described  first  in  the  form  of  an  algorithm,  where  the  elements  described 
above  are  identified  with  steps  in  the  algorithm,  as  are  the  means  to  manage  the  step-wise 
procedure  itself.  Implementation  of  the  procedure  into  a  form  suitable  for  the  production  of 
computational  results  is  described  next.  Example  results  are  presented  for  the  ‘finite  composite 
design’  of  a  cantilevered  beam  under  end  load,  and  for  a  clamped  beam  subject  to  a  distributed 
load.  Sets  of  results  are  obtained  for  both  examples  to  indicate  how  the  layout  of  the  stiffer 
material  (reinforcement)  is  affected  by  varying  the  relative  stiffness  of  the  two  materials.  As 
indicated  earlier,  the  model  for  composite  design  also  is  applied  to  predict  optimal  patterns  for  rib 
reinforcement  of  a  sheet  metal  panel  [  the  interpretation  may  be  contrasted  to  the  treatment  of 
reinforcement  of  plates  in  Diaz  et  al.  (1995)],  and  computational  results  for  an  example  of  this 
application  are  presented  as  well. 


2.  Model  for  the  procedure 

Let  us  consider  the  linear  elastic  structure  occupying  domain  D,  subjected  to  body  forces  f, 
boundary  tractions  t  and  zero  displacement  on  boundary  Fu.  The  structure  will  be  composed  of 
two  materials  identified  respectively  by  their  elasticity  tensors  E1,  for  the  weaker  material  and  E2 
for  the  stiffer  material  (see  Figure  1  a).  In  a  simple  fiber-reinforced  composite,  for  example,  the 
two  moduli  might  correspond  to  the  matrix  and  the  fiber  materials,  respectively. 


The  problem  we  address  here  is  the  following.  Assuming  that  the  domain  Cl  occupied  by  the 
complete  composite  structure  is  specified  and  fixed,  and  given  an  upper  bound  on  the  amount  of 
available  material  '2'  (stiffer  material),  we  seek  to  identify  optimal  properties  and  topolosy  of 
material  s2'  imbedded  in  material '  V  ,  with  no  overlap  of  the  less  stiff  material  in  the  region  of  the 
stiffer  one  (  see  Figure  1  b).  The  objective  for  optimal  design  is  to  maximize  a  measure  of  the 
overall  stiffness  of  the  structure.  The  El  material  properties  are  taken  to  be  specified  and  fixed, 
with  E’>  0.  With  Cl2  to  symbolize  the  part  of  the  domain  taken  up  by  concentrated  material  ‘2’  (to 
be  optimally  designed),  it  is  required  that  in  the  final  design  -  Clz.  Since  the  structural 

domain  Cl  is  fixed,  this  indicates  that  material  ‘Y  fills  the  part  of  the  original  domain  not  occupied 
by  the  optimal  material  ' 2 '  without  overlap. 


a)  Initially  as  a  perfect  mixture.  b)  As  a  two  material  composite. 

Figure  1-  The  structure  subjected  to  body  force  f  and  boundary  traction  t. 

In  order  to  solve  the  problem  described  above,  the  following  procedure  is  proposed: 

Initial  design 

As  a  first  step  assume  that  one  has  a  perfect  mixture  of  the  two  materials  (i.e.  the  effective 
material  tensor  of  the  mixture  is  given  by  Eeff=E'+E2;  note  that  for  the  2D  model  to  simulate 
laminated  structure  [see  e.g.  Pedersen  (1993)]  the  'mixture'  provides  an  authentic  model  for  the 
effective  material  modulus).  Again,  recall  that  at  the  final  design  the  materials  are  to  be  effectively 
separated. 

For  this  mixture  and  identifying  the  trace  of  the  E2  tensor  as  the  sensible  measure  of  the 
material  (  p  =  rr(E2)),  design  material  to  obtain  its  optimal  local  properties  and  distribution  in 
Cl.  Here  the  resulting  design,  obtained  using  the  method  described  in  Bendsoe  et  al.  (1994),  is 


represented  by  continuously  varying  p  (  the  plot  of  p  is  sometimes  referred  to  as  the  'shades  of 
gray'  diagram).  Once  this  solution  for  the  optimal  material  and  its  distribution  over  the  structure 
are  known,  a  finite  number  of  optimization  sub-problems  are  performed  to  predict  a  refinement  of 
it  into  a  design  having  total  material  separation  and  concentration  of  material  '2'  at  its  upper 
bound.  The  redesign  of  material  '2'  is  accomplished  using  a  method  (Guedes  and  Taylor  1997  a) 
for  the  prediction  of  optimal  topology.  The  method  is  applied  stepwise,  where  in  each  step  a 
gradually  higher  value  of  unit  cost  is  ascribed  to  regions  of  relatively  low  value  of  p.  Specific 
details  are  provided  next,  where  the  procedure  is  first  described  formally  in  the  form  of  an 
algorithm,  and  is  later  expressed  for  computational  treatment  using  a  finite  element  interpretation 
of  the  continuum  structure. 


2.1.  Algorithm 

It  is  assumed  that  an  initial  design  (the  shades  of  gray  results  described  above)  has  been 
obtained.  Based  on  this  result  and  recalling  that  ’design’  is  represented  by  the  p  field,  the  algorithm 
is  described  as  follows: 


Step  1:  Define  an  increasing  sequence  of  N  evenly  spaced  cutoff  values  pek ,  k=],2...N  such 
c  p  _  c  Tp') 

that  p[  =  —  and  pcN  =  p  ;  pct  =  k  —  .An  optimization  sub-problem  is  to  be  solved  for  each 
N  l  N  J 

value  of  pck ,  where  index  k  identifies  steps  in  the  procedure.  The  number  N  designates  the  number 
of  steps  to  achieve  the  final  design.  Both  N  and  the  upper  bound  p  are  prescribed.  Considerations 
of  how  properly  to  select  values  of  these  data  are  discussed  below. 


Step  2:  For  a  given  cut  off  value  pck ,  identify  the  following  sub-domains  of  the  structure: 


-  {x  €  D. :  pk  >  pckj,  sub-domain  with  higher  value  of  material  measure  p  ( 1 ) 

D"  =  {xe  Q  :pk  <  ptc},  sub-domain  with  lower  value  of  material  measure  p  (2) 


Qt  =  {xe  O  :pk  =  p},  sub-domain  with  material  measure  p  at  its  upper  bound.  (3) 


Recall  that  p  is  identified  as  the  trace  measure  of  the  elasticity  tensor  E2 .  and  note  that 
Clk  c Qk  and  Q*  uQ~  =  Q. 


Step  3:  With  the  objective  of  driving  the  designate  material  at  each  point  to  either  its  upper 
(p)  or  lower  (p)  bound  value,  a  unit  relative  cost  coefficient  <uis  modified  stepwise  in  a  manner 

to  exaggerate  the  cost  of  material  wherever  p  has  value  below  the  cutoff.  Specifically,  the  relative 
unit  cost  function  based  on  the  results  at  step  k  is  set  equal  to  the  constant  value  one  or 
5)  according  to 


[  1  for  x  e  D.*k  ' " 

<4> 

with  value  of  the  constant  co  »  1 . 

Total  cost  at  the  (k  +  1  )s[  step  is  evaluated  according  to: 

\cokPk^dD.  .  (5) 

0 


Step  4:  Once  the  relative  unit  cost  is  defined,  optimize  the  material  ‘2’  elasticity  tensor  E2 , 
i.e.,  determine  the  measure  pi4.,  and  the  remaining  attributes  of  E2  using  the  cost  distribution 

determined  in  step  3.  This  is  accomplished  by  solving  the  (shades  of  gray  )  problem: 

E,  Min  (6) 

o<psp,.,»rr(E:)sp  *-n  r  > 

subjected  to  the  following  constraints: 

J  (l-*(oj£;,  ^(u)  ea(v)rfn  +  j£^  etj (u)  ekl(v)dCl- 

n  q 

-j/,v,dn-J(,virfr=o 

n  r, 


Vvs  V 


(7) 


(S) 


j  pwrfn  +  j  mp^dn<  r. 

n:  n; 


where  V  represents  the  set  of  admissible  displacements  and  where  x  (p  t  )  is  the  characteristic 
function  defined  in  terms  of  the  set  H  4  as: 


*(&)« 


f  1  for  x  e  D ,  1 
{ 0  for  flj 


(9) 


Note  that  u  and  E2  in  (7)  are  intended  to  stand  for  the  solution  of  the  current  step. 
An  alternative,  equivalent  statement  of  the  equilibrium  condition  (7)  is: 


Jfe/  +£i )«*(“)  etl(v)dQ  +  j E?jtl  ei;.(u)  eu(v)^n- 


n-nt  nt 

-  J  /,  v,.  -  Jr;  v;  *£T  =  0  V  v  e  V 

n  r, 


(10) 


Note  that  in  form  (10)  for  the  problem,  separation  of  the  structure  in  Q  into  the  domain 
Clt  where  p(x)  =  p  and  the  remaining  Q-Qt  is  reflected  explicitly  in  the  integrals. 


The  lower  bound  p  is  set  to  have  the  value  on  the  order  10‘7  p  . 

For  the  equilibrium  constraint  (7)  note  that  the  fixed  material  T  exists  (  mixed  with  material 
'2')  only  where  the  measure  of  material  '2'  is  not  its  upper  bound  i.e.  material  T  is  removed  from 
the  region  with  full  material  '2'.  As  the  procedure  converges  this  will  provide  full  pointwise 
separation  of  the  two  materials.  Also  in  the  resource  constraint  we  assume,  as  previously  stated, 
that  the  material  unit  cost  for  higher  value  of  the  material  measure  p  is  equal  to  one,  and  for  the 
lower  value  of  the  measure  is  to  »1.  The  distribution  of  material  T  is  fixed  during  this  step,  and 
the  solution  obtained  provides  the  description  pkM  of  the  trace  of  (E2)k+i  as  well  as  the  current 

local  material  properties  (e,2^  )  . 


Step  5:  Set  the  cut  off  value  to  p*+1 

Go  to  step  2  and  repeat  the  procedure  until  the  limiting  cutoff  value  p  is  reached. 


2.2.  Optimal  distribution  of  the  designable  material 

We  will  analyze  in  more  detail  the  optimization  subproblem  to  be  solved  in  step  4  of  the 
previous  section.  The  aim  of  this  step  is  to  solve  the  following  problem:  For  a  prescribed  material 
unit  cost  function  <y(x)  and  fixed  sets  Cl..,  Cl[  andf^,  find  the  optimal  distribution  of  designable 
material  as  the  one  that  minimizes  the  structural  compliance  (maximizes  the  stiffness)  of  the 
structure,  as  stated  in  the  optimization  problem  (6-8). 

Following  the  presentation  in  the  cited  Bendsoe  et  al.  (1994)  paper  for  the  treatment  of  this 
problem,  the  optimization  problem  can  be  equivalently  stated  as, 


subjected  to  the  resource  constraint , 

jtr(E2)dH  +  j5Jtr(E2)dD<R.  (12) 

n;  n- 

Note  that  in  this  statement  of  the  structural  compliance  optimization  problem,  the  equilibrium 
constraint  is  interpreted  in  the  form  of  a  minimum  potential  energy  statement. 

The  design  variable  E2  is  characterized  by  the  field  p(x)  which  identifies  global  distribution  of 
the  trace  measure  of  E2,  and  by  its  local  tensorial  structure.  The  maximization  problem  in  (1 1)  can 
be  split  into  two  successive  maximization  problems  (see  Bendsoe  et  al.  1994),  associated 
respectively  with  these  'field'  and  'local  structure'  attributes  of  E2. 

Thus  (1 1)  is  rewritten  as: 


where  the  inner  maximization  identifies  the  optimal  elasticity  tensor  within  a  prescribed  value  of 
the  trace  p  (optimal  relative  values  among  tensor  components)  and  the  outer  max  establishes  the 


optimal  distribution  of  the  trace  of  the  elasticity  tensor  within  the  prescribed  upper  and  lower 
bounds  (optimal  material  distribution). 


From  the  fact  that  the  resource  constraint  is  only  dependent  on  the  design  variable  p  ,  /. e. ,  on  a 
norm  of  the  tensor  E  (the  trace  in  our  case)  and  not  on  its  individual  components,  consider,  for 
now.  only  the  inner  max  min  problem  in  (13).  Following  the  same  argument  as  presented  in 
Bends0e  et  al.  (1994),  the  inner  max  and  min  can  be  interchanged  as, 


Min  Max 

v€V  £*») 

"(e2)=p 


«*(v)</n+f£*,«  (v)  eu(v)dQ. 

n  n 

“  J  f;vi  dCl-j t,v,  dT 


(14) 


Following  this  transformation,  the  maximum  can  be  solved  analytically  for  the  design  variable 
E:  and  has  the  solution,  ‘  ‘  • 


rv  _  «,y(v)eu(v) 

Eijkl  “  P  * 


Mv>f 


with  the  optimal  strain  energy  density  equal  to, 
ea{v)=  ipe{(v)  e#(v)  . 


Substituting  (16)  into  problem  (13)  one  obtains, 


Max  Min  j 

0  SpSp  v€V 


Jt1  -  x(f^k  ))h Jae-(v)  eu  (v)d £l  +  \p  5ik5j(  e;j ( v)  ekl  (v) dft 
n  n 

-Jf.v.dn-JtjVjdr 


subjected  to  the  resource  constraint , 

J  pdQ  +  jciJ  pdD.  <  R  . 

Ot  n£ 


(15) 


(16) 


(17) 


(IB) 


(Note:  Result  (15)  is  valid  only  for  one  load  condition). 

The  previous  problem  is  a  maxmin  problem  where  the  inner  minimization  has  a  unique 
solution.  Thus  the  problem  is  differentiable  and  the  necessary  conditions  for  a  maximum  are, 


where  u(x)  is  the  unique  solution  of  the  inner  minimization  i.e.  the  displacement  of  the  structure  at 
equilibrium.  Note  that  these  conditions  are  to  be  satisfied  at  the  optimum. 


In  the  previous  conditions  A  is  the  Lagrange  multiplier  for  the  resource  constraint  and  r\u  and 
r\i  are  the  multipliers  for  the  upper  and  lower  bound  constraints  respectively. 

Complementarity  conditions  for  the  problem  are: 

T7„(p-p)  =  0,  TJU  >0,  p  <p  VxcH  (21) 

77,(p-p)  =  0, 77,  >0,  p  >p  Vxsfl  (22) 


A 


J  pdD.  +  J  co  pdO.  -  R 

n;  n; 


=  0,  A  >  0 


(23) 


It  follows  from  (21,  22)  that 

tiu71,=0  Vxe  H  (24) 

From  the  preceding  conditions  it  can  be  shown  that  A  >  0,  implying  that  the  resource 
constraint  will  be  active. 


2.3.  Optimal  distribution  of  the  designable  material  (alternative  formulation) 

As  an  alternative  to  the  part  of  the  procedure  of  section  2.1  where  the  fixed  material  is 
removed,  it  is  possible  to  reformulate  the  problem  of  optimal  distribution  of  material  so  that  the 
fixed  material  is  gradually  removed  within  the  process  of  optimizing  the  material  distribution. 

With  this  purpose  in  mind,  it  is  assumed  that  within  the  material  distribution  optimization,  the 
domain  occupied  by  the  fixed  material  T  is  the  complement  of  the  domain  H(p)=  {x  :p  =  p}, 

and  consequently  a  function  of  p.  Thus  problem  (13)  is  restated  has, 


\  J  t1  -  Xt  (p(p)))Eweij  (' V)  (v)rfn  +  j e..  (v)  eu  (v)dQ. 

“U  n 

-jfivid&-jtividr 

n  r, 

where  the  *t(Q(p))  is  an  £  differentiable  regularization  of  the  characteristic  function  of  the  set 
&{p)  =  {x  e  12 :  p  =  p},  having  the  properties  (  see  Figure  2): 

if  p  =  p 

if(l-£)p<p<p  (26) 

if  p  <  (l-e)p 


1 


Figure  2-  Regularization  of  the  characteristic  function. 

Note  that  this  regularization  is  non-symmetric  in  the  e  neighborhood  of  p  and  thus  imposes  the 
effective  separation  of  the  two  materials  on  the  region  £2(p)  (/.e.,  at  all  the  points  where  p  —  p). 

Following  the  rationale  of  the  previous  section  we  will  have  the  same  characterization  of  the 
optimal  material  E  (15),  but  the  necessary  conditions  on  p  (  written  for  the  entire  domain)  are, 

EUea  (u)ew  (u)  +  ei}  (u)e?(u)  -  A -T)u  +p,  =0  Vxe  Q. 


^xAp) 

dp 


1 

,  v  _  (p  +  £p-p):(gp+2p-2p) 

‘  (pej* 

0 


Max  Max  Min 

0SpSpE->0.  vsV 
tr{z-)=p 


P  =  P 


(27) 


where,  A,  r\u  and  77 ,  satisfy  the  orthogonality  conditions  (21-24)  and 


^Xc(p)  _ 

dp 


6(p-p  +  £p)(p-p) 

fa? 


if  p  =  p 


if  (l-£)p  <  p  <  p 


if  p<(l-g)p 


(28) 


The  purpose  of  this  a  regularization  is  to  implicitly  remove  the  fixed  material  E1  in  regions 
within  D.  where  the  designate  material  is  at  its  upper  bound.  Note  that  for  every  £  ?;  0  the 
problem  as  stated  in  (25)  is  a  meaningful  design  problem.  However,  the  limit  case  (/.«.,  when  £=0) 
is  not  a  well  posed  problem,  in  the  sense  that  one  can  get  as  close  as  one  wants  to  the  designable 
material  upper  bound  and  still  have  fixed  material  coexisting  with  the  designable  material. 


3.  Method  for  computational  solution 

As  described  above,  the  final  structural  topology  is  predicted  through  an  iterative  procedure 
where  the  designable  material  is  concentrated  at  its  upper  bound  p  (remember  that  the  trace  of 
elastic  tensor  is  identified  as  such  a  measure).  This  is  attained  through  the  solution  of  a  series  of 
material  optimization  sub-problems  (defined  in  section  2.2)  where  the  material  unit  cost  is 
iteratively  adjusted  so  that  lower  trace  measures  impose  a  higher  cost  on  the  resource  distribution. 

In  the  following,  the  implementation  of  these  steps  for  computation  is  discussed. 


3.1.  Material  unit  cost 

The  basic  idea  behind  a  non-uniform  material  unit  cost  imposed  on  the  resource  constraint  is  to 
induce  a  concentration  of  the  designable  material  at  its  upper  bound  (full  material).  As  the  first 
step,  the  optimization  problem  is  solved  initially  assuming  a  uniform  material  unit  cost.  Based  on 
results  for  this  optimal  solution,  the  structural  domain  is  divided  into  two  regions,  the  part 
n  (below  cutoff  value  pc)  and  D.  (the  part  above  p‘),  as  defined  earlier.  Once  these  sets  are 


defined,  one  computes  the  new  optimal  material  distribution  imposing  a  higher  unit  cost  on 
material  in  Cl~ ,  i.e.  the  resource  constraint  is  now  changed  to, 


|  p  dCl  +  |  co  p  dQ.  <  R,  co  »  1 

cr  cr 


(29) 


In  figure  2  this  process  is  simulated  in  a  one-dimensional  sketch  for  initial,  intermediate  and 
final  distributions.  Assuming  that  the  material  distribution  is  known  for  the  initial  design,  the 
resulting  distribution  p(x)for  the  subsequent  steps  (general  step  k  and  final  step  N)  reflects  a 

relative  increase  value  in  regions  (fl -Q").  The  increase  for  the  k‘hstep,  and  the  concentration  of 
material  at  the  upper  bound  for  the  final  design  are  shown  in  separate  curves. 


Figure  3-  Sketch  of  material  measure  p(x)  for  initial,  intermediate  and  final  designs  (typical). 

The  idea  behind  this  approach  is  simple,  but  one  major  feature  of  the  method  needs  to  be 
addressed.  It  follows  from  the  fact  that  the  procedure  attributes  higher  unit  costs  according  to 
regions  occupied  by  low  material  measure  p,  rather  than  to  the  material  itself.  As  a  consequence, 

once  a  region  is  included  in  the  Cl  sub-domain  it  tends  to  remain  there.  To  deal  with  such  a 
feature,  care  should  be  taken  in  the  selection  of  the  number  of  cutoff  values. 


3.2.  Optimal  material  distribution 


Once  the  material  unit  cost  function  (&)k  at  the  k'h  step  of  the  topology  design  procedure)  is 
defined,  one  can  calculate  the  new  optimal  material  distribution.  Note  that  the  domain  occupied  by 
the  fixed  material  is  known  and  fixed  during  this  step. 


The  methodology  adopted  is  based  on  a  sequential  solution  of  the  optimality  conditions  (19, 
20)  with  the  strain  field  calculated  using  a  finite  element  approximation  of  the  elastostatic  problem, 
characterized  by  the  equilibrium  statement  (7)  (the  commercial  program  ANSYS  was  exploited  for 
this  purpose). 


So  based  on  the  material  unit  cost  for  the  procedure  step  k  (<Dk)  and  assuming  the  design 
variable  p'  constant  at  each  finite  element,  one  can  write  the  optimality  conditions  separately  for 
each  element.  From  the  discrete  interpretation  of  the  these  conditions,  and  introducing  the  upper 
and  lower  bound  constraint  thickness  parameter  C,  (defined  by  the  user),  the  solution  is  obtained  by 
the  fixed  point  method. 


max  fo-OfcOi-p]  if  a'fcOj  <  max  [0-cXp*i*p] 

(p'l,  =<  a'(p'\  ifmaxfo-cXpOi.pka'fcOj  £  min[(l  +  £)(pe).,p], 
™n[(l  +  £)(p*)..p]  if  min[(l  +  £)(pe).,p]<ae(pe). 


(30) 


with  the  multiplier  ae  given  by, 


a 


_  (eijfrW), 

2co,A 


for  each  finite  element. 


(31) 


In  the  previous  algorithm,  the  index  'e'  ranges  over  all  the  finite  elements,  T  is  the  iteration 

counter  and  <  >e  identifies  the  average  operator  applied  to  the  eth  element  strain  field. 

The  calculation  of  the  Lagrange  multipliers  A  is  determined,  within  each  iteration,  through  a 
updating  scheme  that  imposes  the  resource  constraint. 


=  R. 


(32) 


\ 

j  cop  dH-r  j  pdD 


In  the  case  of  the  alternative  formulation  described  in  section  2.3,  where  the  domain  occupied 
with  fixed  material  is  gradually  removed,  the  fixed  point  algorithm  is  as  previously  stated  in  (30) 
with  the  of  factor  given  by. 


cr‘  = 


(eg(u>B(u))e 


tylcyi  +  ^-EikI(eij(u)eij(u))4 


(33) 


4.  Examples 

The  examples  presented  below,  try  to  demonstrate  the  applicability  of  the  developments 
described  within  the  context  of  three-dimensional  structural  applications.  The  first  two  examples 
are  discretized,  for  computational  purpose  with  eight  node  isoparametric  brick  elements,  an 
uniform  finite  element  mesh  and  the  computational  prbeedure  is  programmed  within  the  Ansys 
finite  element  software. 


4.1.  Example  1 

The  first  example  considered  is  the  problem  of  finding  the  optimal  topology  of  material  ‘2’  in  an 
end  loaded  cantilever  beam  made  of  an  isotropic  material  with  Young  Modulus  El  and  Poisson 
coefficient  0.2  (material  ‘1’).  The  design  domain  has  dimensions  20x12x1  (see  Figure  4)  and  the 
resource  constraint  upper  bound  equals  40%  of  the  volume  of  the  design  domain. 

Material  two  has  Young  modulus  upper  bound  of  109  Pa  and  a  lower  bound  of  101  Pa.  The 
black  and  white  procedure  takes  N=40  steps  and  the  shades  of  gray  optimization  (step  4  of  the 
procedure)  is  limited  to  5  iterations. 

This  example  is  the  3D  equivalent  of  the  2D  example  described  in  Guedes  and  Taylor  (1997). 
Even  though  three  dimensional,  the  model  would  behave  like  a  two-dimensional  one  due  to  the 
small  design  domain  thickness.  This  fact  will  permit  the  comparison  of  the  results  obtained  with 


the  ones  presented  for  2D  solutions  in  the  case  of  non-existent  material  *1’  (see  e.g.  Guedes  and 
Taylor  (1997)). 

The  solutions  obtained  are  shown  for  various  ratios  between  the  maximum  Young  modulus  for 
material  ‘2’  (max  E2,  attained  at  p  =  p )  and  the  Young  modulus  for  the  fixed  material,  E1. 

The  first  solution  has  a  small  value  of  this  ratio  (10’7)  and  consequently  is  equivalent  to  the 
optimal  topology  design  problem  restricted  to  one  material.  The  subsequent  results  consider 
higher  material  ratios. 


Figure  4-  Design  domain  and  boundary  conditions. 


a)  b) 

Figure  5-  Optimal  topology  for  E'/(max  E2)=10'7. 


a) 


b) 


Figure  6-  Optimal  topology  for  El/(max  E2)=l/50. 


a)  b) 

Figure  7*  Optimal  topology  for  E‘/(max  E2)=  1/10. 


Figures  on  the  left  show  the  designable  material  imbedded  within  the  fixed  material.  In  the 
right,  and  to  better  portray  the  final  topology  of  the  designable  material,  the  fixed  material  is  not 
shown.  Note  also  that  for  representation  purposes  the  designable  material  (material  *2’)  is 


identified  in  black  color  when  the  two  materials  are  present  and  in  lighter  gray  when  material  ‘T  is 
removed  from  the  figure. 

From  the  results  obtained  one  can  observe  that  the  truss  like  structure  inside  the  domain 
weakens  as  the  material  ratio  increases.  This  is  attributed  to  the  fact  that  as  material  T  is 
stiffened,  it  has  an  higher  contribution  in  carrying  the  shear  load  present,  thus  freeing  material  ‘2’ 
to  support  the  bending  load. 


4.2.  Example  2 

This  example  considers  the  topology  optimization  of  material  ‘2’  for  the  cantilever  beam 
presented  in  Figure  9.  The  design  domain  has  dimensions  8x8x20  and  is  discretized  with  10264 
eight  node  isoparametric  brick  elements.  The  applied  force  is  restricted  to  a  shear  force  distributed 
as  shown  below.  The  resource  constraint  upper  bound  is  eQual  to  40%  of  the  design  -domain 
volume.  Material  ‘2’  has  a  Young  modulus  upper  bound  of  109  and  a  lower  bound  of  101.  The 
black  and  white  procedure  was  done  in  N=20  steps  and  the  shades  of  gray  optimization  (step  4  of 
the  procedure)  is  limited  to  5  iterations. 


Figure  9  -  Design  domain  and  boundary  conditions. 


The  figures  below  display  the  results  obtained  for  different  ratios  of  E‘/(max  E:).  For  this 
example,  the  fixed  material  (  ‘1’)  is  assumed  isotropic.  For  each  material  ratio,  the  final  result 
internal  structure  can  be  visualized  in  the  various  transversal  and  longitudinal  cuts  shown. 


Figure  10  -  Final  topology  for  E'/(max  E2)=10‘\ 
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Figure  13-  Final  topology  for  E'/(max  E‘)=  1/10  (longitudinal  cuts). 


Figure  14-  Final  topology  for  E'/(max  E2)=  1/3. 
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Figure  15-  Final  topology  for  E'/(max  E:)=l/3  (longitudinal  cuts). 


For  this  last  material  ratio,  the  next  figure  shows  the  final  topology  of  material  *2’  with 
material  ‘  1  ’  removed. 


Figure  16-  Final  topology  for  E'/(max  E2)=l/3  (material  one  not  shown). 


Figure  17-  Computer  model  of  the  car  decklid. 


Figure  18-  Optimum  topology  of  embossed  ribs  using  30%  of  the 


Note:  For  representation  purposes  the  designable  material  (material  ‘2’)  is  identified  in  black  color 
when  the  two  materials  are  shown  and  in  lighter  gray  when  material  T  is  removed  from  the 
figure. 

We  note  here  the  big  gap  between  the  first  (10'7)  and  second  material  ratio  (10‘‘)  used  in  the 
examples.  This  is  due  to  the  fact  that  the  changes  in  topology  as  function  of  the  material  ratio  are 
only  easily  detectable  for  ratios  above  10’\ 


4.3.  Example  3  Embossed  Ribs  in  Stamped  Plates 

The  design  model  described  in  this  work  has  many  applications,  among  them,  the  optimum 
design  of  reinforcing  embossed  ribs  (also  known  as  beads)  in  stamped  plates.  A  particular  case  of 
the  optimization  problem  (13)  is  obtained  when  the  local  structure  of  the  material  '2'  is  fixed  and 
orthotropic,  and  when  material  T  is  isotropic.  In  such  cases  it  is  possible  to  determine  the 
optimum  layout  of  embossed  ribs  (material  7")  within  the  flat  plate  (material  T).  In  order  to  solve 
this  particular  optimization  problem  the  only  modification  required  is  to  fix  the  form  of  E2  in  the 
most  inner  maximization  problem  in  (13). 

Due  to  the  geometric  complexity  of  this  real  life  example  it  was  solved  using  software 
developed  at  Ford  for  this  particular  application  of  optimum  layout  of  embossed  ribs  of  constant 
properties.  The  structure  in  question  is  the  decklid  of  a  sedan  vehicle  (see  Figure  17)  subject  to  a 
vertical  load  at  the  key  switch.  The  objective  of  the  problem  is  to  maximize  the  stiffness  for  such 
load  condition  by  embossing  ribs  in  the  decklid  inner  panel.  The  resource  constraint  metric  is  the 
area  allowed  to  be  covered  by  embossed  ribs.  For  this  example,  embossed  areas  are  limited  to  30% 
of  the  total  decklid  area.  Figures  18  and  19  show  the  ‘optimum  topology’  and  the  ‘optimum 
orientation  of  the  principal  directions’  of  the  embossed  ribs,  respectively. 

The  stiffness  of  the  decklid  was  increased  100%  with  the  resulting  topology/orientation  of  ribs. 


5.  Discussion 


To  contrast  the  present  treatment  of  two-component  composite  design  with  common  models 
for  the  design  of  continuum  structures  based  on  homogenization  interpretation  of  a  two-phase 
materials  [the  subject  is  surveyed  in  Bendsde  (1995);  see  also  Cherkaev  and  Gibiansky  (1997); 
Lipton  (1988)],  note  that  in  the  latter  the  local  properties  are  designated  to  have  specific  form,  e.g. 
isotropic,  rather  than  to  be  represented  by  a  free  modulus  tensor.  Also,  the  'zero-one' 
interpretation  for  the  prediction  of  optimal  topology  from  the  latter  models  [Rozvany  (1995)  et  al. 
provides  a  survey  on  methods  for  topology  design]  results  in  a  ill-posed  problem,  whereas  in  the 
present  approach  topology  is  predicted  on  the  basis  of  an  ordered  sequence  of  solutions  to  a  well- 
posed  optimization  problem. 

For  simplicity,  the  model  for  two-component  composite  design  is  described  in  this  paper  in  a 
relatively  narrow  form.  However,  the  same  method  may  be  applied  equally  well  in  other  contexts- 
for  structural  design  with  non-linear  materials  as  represented  in  Bendsde  et  al.  (1996)]  for 
example,  or  for  problems  where  the  local  measure  of  material  cost  is  expressed  by  an  invariant 
other  than  the  trace  of  the  modulus  tensor  [see  e.g.  Taylor  and  Washabaugh  (1995),  Taylor  (1998) 
for  expressions  of  generalized  cost]. 
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